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^ Abstract 

Some of the most impressive singular wave fronts seen in Nature are 
the transbasin oceanic internal waves, which may be observed from 
I the Space Shuttle as they propagate and interact with each other, 

^ for example, in the South China Sea. The characteristic feature of 

these strongly nonlinear wavefronts is that they reconnect when two 
of them collide transversely. We derive the EPDiff equation, and use it 
to model this phenomenon as elastic collisions between singular wave 
fronts (solitons) whose momentum is distributed along curves moving 
in the plane. Numerical methods for EPDiff based on compatible dif- 
ferencing algorithms (CDAs) are used for simulating these collisions 
among curves. The numerical results show the same nonlinear behav- 
ior of wavefront reconnections as that observed for internal waves in 
the South China Sea. We generalize the singular solutions of EPDiff 
for other applications, in computational anatomy and in imaging sci- 
ence, where the singular wavefronts are evolving image outlines, whose 
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momentum may be distributed on surfaces moving though space in 
three dimensions. The key idea is always momentum exchange during 
coUisions of the wavefronts. A suite of 2d and 3d numerical simulations 
provide collision rules for the wavefront reconnection phenomenon in 
a variety of scenarios. 
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1 Synopsis 

Space Shuttle observations show great lines in the sea, which show up in Syn- 
thetic Aperture Radar (SAR) images as in Figure[l| These are internal waves, 
whose fronts appear in SAR images as unbroken curves extending for hun- 
dreds of kilometers. For example, tidal flows and undulations of the Japanese 
Current in the region of the Luzon Strait between Taiwan and the Philippines 
generate internal wavefronts which are well over one hundred kilometers in 
length and which propagate westward for hundreds of kilometers all the way 
across the South China Sea. These wave fronts show strong coherence and 
possess the characteristic feature of reconnecting with each other whenever 
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any two of them intersect transversely. This reconnection is the hallmark 
of a nonlinear process. Weaker waves may intersect, but not reconnect. We 
seek to encode the motion of these internal wave fronts and their collision 
rules mathematically in a minimal PDE (partial differential equation) model, 
and to investigate these wave front reconnections in numerical simulations 
of the model. Establishing a simplified PDE model whose solutions encode 
the motion of these internal wave fronts and their collision rules will enable 
understanding the effects of the nonhydrostatic processes that govern them, 
without requiring the full numerical simulation and analysis of the 3d fluid 
motion equations. 

We shall also derive more detailed PDE models for the interactions of 
these nonlinear wavefronts when topography and boundaries are also in- 
cluded. However, we shall defer developing the numerical methods needed 
to deal with these additional effects until elsewhere. 

Our derivation of the 2d multilayer Euler-Poincare (EP) equations be- 
gins by vertically integrating the variational principle for Euler's equations 
in 3d. These 2d equations must be nonhydrostatic, because the wave fronts 
they model possess strong horizontal gradients of vertical acceleration. These 
nonhydrostatic terms are included in a hierarchy of EP equations which is 
derived by making a series of simplifying approximations in the exact vari- 
ational principle. The resulting multilayer 2d equations preserve the EP 
properties of energy balance, circulation laws and potential vorticity (PV) 
conservation. PV analysis in a multilayer fluid system is essential for as- 
sessing its baroclinic instability [59]. However, PV for wave fronts is a new 
concept, whose implications are hardly understood yet. We will show that 
the new hierarchy of EP equations recovers a sequence of previously known 
models of internal wave dynamics when specialized in various ways. 

Next, from the full hierarchy we extract a simple, minimal PDE descrip- 
tion, called EPDiff, which models internal wave fronts as delta functions of 
momentum distributed on moving curves in the plane. This corresponds to 
modeling internal wave fronts as contact discontinuities in velocity. EPDiff 
thus explains the creation and stability of wave fronts as the development 
of singular momentum solutions from continuous velocity distributions in a 
PDE initial value problem (IVP). This description links the shape of the 
wavefront velocity profiles to the Green's function associated with the non- 
local relation in the PDE between its singular momentum and its continuous 
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velocity. Hence, EPDiff explains the uniform width of the wave front velocity 
profiles that is observed in internal wave trains. Namely, the momentum of 
the wavefront is concentrated on delta functions supported on the curved 
evolving wavefronts. The corresponding velocity profile is then obtained via 
a Green's function relation between wave momentum and fluid velocity. This 
minimal description of internal wave fronts as singular momentum solutions 
of the evolutionary equation EPDiff encompasses their reconnections and 
provides the geometric mechanism underlying their propagation and colli- 
sion interactions. In Id, these collisions are understood as elastic solition 
collisions, whose solutions are obtained by the inverse scattering transform 
(1ST) for an associated isospectral linear eigenvalue problem. As far as we 
know, the machinery of 1ST is not available in higher spatial dimensions for 
EPDiff. However, the elastic scattering interactions seen in its wavefront 
solutions in 2d (and in 3d) are explained by another geometric mechanism. 
As we shall explain, this geometric mechanism for propagation and collisions 
turns out to be Hamiltonian geodesic flow on the time-dependent smooth 
maps (diffeomorphisms), defined with respect to the kinetic energy metric 
that also determines the wave profile. 

To solve the EPDiff equation, we develop a numerical method called the 
compatible differencing algorithm (CDA) that is able to capture the collisions 
among these weak solutions of EPDiff, and we characterize the wave front 
interactions we observed numerically in a variety of scenarios. 

In summary, we use EP variational theory to derive the EPDiff equation, 
whose solution shows singular wave fronts, and we use new CDAs for its 
numerical simulations. EPDiff is the first mathematical explanation of the 
observed 2d internal wave front reconnections. 

Geometric setting. We have found that modeling the reconnection pro- 
cess in internal wave front dynamics requires a class of PDEs that are both 
nonlinear and nonlocal. EPDiff is a not exactly a hyperbolic equation. It has 
a characteristic velocity, but the relation its EP variational principle defines 
between the fiuid's velocity and momentum in Newton's Law is nonlocal. 
That is, the fluid's velocity is determined from its momentum by solving 
an elliptic equation. Physically, the elliptic equation arises from nonhydro- 
static processes that cause linear, or nonlinear, dispersion or focusing. This 
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nonlocal relationship between velocity and momentum is reminiscent of the 
Biot-Savart relationship between velocity and vorticity. However, like soli- 
tons, internal wave fronts carry momentum and inertia, while fluid vortices 
do not. 

Our work has some close similarities with soliton theory and some im- 
portant differences from it. The Korteweg-de Vries (KdV) equation at linear 
order in the asymptotic expansion for shallow water waves and the Camassa- 
Holm (CH) equation at quadratic order are both soliton equations, and they 
are both associated with geodesic flow [33l [SU HT] . That is, they each make 
optimal use of their kinetic energy, which provides a norm for their velocity. 
The EPDiff equation has this same characteristic. Moreover, the singular so- 
lution ansatz for the geodesic flow of momentum associated with the EPDiff 
equation which we discuss below turns out to be a momentum map, as dis- 
covered in Holm and Marsden [22] . The momentum map property of these 
singular solutions means they comprise an invariant manifold, preserved by 
the flow of EPDiff. This property allows us to reduce the dimension of the 
wave front interactions to an invariant manifold of singular momentum so- 
lutions of EPDiff. These solutions describe the observed propagation and 
reconnection phenomena of the wave fronts, but they have no internal de- 
grees of freedom, and thus they have no mechanism for wave breaking to 
occur. 

In this paper, we model internal wave fronts as contact discontinuities 
in velocity, whose motion is governed by applying time-dependent smooth 
maps which act on the curves in the plane that outline the wave fronts and 
specify their momentum vectors at points along these curves. The motion of 
the internal wave fronts is governed by the EPDiff equation. This equation 
is the condition for the smooth maps that act on these curves to evolve by 
geodesic flow, with respect to a metric determined by their kinetic energy, 
which also determines the functional form of their wave profile. The wave 
front motion may thus be determined numerically as an initial value problem 
by solving the EPDiff equation. 

In such a strongly geometric setting, we expect that combining finite 
element methods and discrete exterior calculus may provide future improve- 
ments in both modeling and simulations, by providing a useful setting for 
integration of the analysis with the numerics in the description of internal 
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wave fronts as singular weak solutions, or limiting solutions (contacts) for 
nonlinear, nonlocal PDEs. However, we defer the full use of discrete exterior 
calculus methods for numerically solving EPDiff until these methods have 
been further developed. See [19], [35] for introductions to these methods. 



Outline of the paper. Mathematical modeling problems for internal wave 
fronts, and our approaches to solving them, are summarized in Section |2} A 
preview of our numerical results appears in Section [3j Previous work in mul- 
tilayer descriptions of internal waves is reviewed in Section |4j The columnar 
motion ansatz for multilayer fluids is explained in Section |5j and is used there 
to derive the nonhydrostatic 2d MultiLayer Columnar Motion (MLCM) equa- 
tions and their natural boundary conditions. A series of weakly nonlinear 
approximations is introduced into the EP variational principle for the MLCM 
equations in Section |6} The resulting weakly nonlinear 2d multilayer equa- 
tions are also compared to the Id multilayer equations of Choi and Camassa 
pT] and others. In Section [7| we derive the EPDiff geodesic equation by 
neglecting dispersion due to potential energy. Singular solutions of EPDiff 
and their canonical dynamics in the Lagrangian fluid representation are also 
introduced in Section [T] 

The main numerical results of this paper are given in Sections [8} [9] and 
10 Section [8] describes our numerical approach using compatible differencing 



algorithms. Section [9] provides a large suite of numerical simulations investi- 
gating the collision rules for the interactions of wave-front solutions of EPDiff 
in 2d. Section 10 extends our numerical solutions of EPDiff to 3d, and shows 
that its stable codimension-one singular solution behavior persists in higher 
dimensions. Finally, we discuss conclusions, future directions and remaining 



outstanding problems for internal wavefront interactions in Section 11 



2 Problem statement, approach, and main 
results 

Space Shuttle observations of the South China Sea surface show a sequence 
of large amplitude internal waves at basin scale moving westward after being 
created by tides and currents running through the Luzon Strait on the eastern 
side of the basin. These internal waves appear in Synthetic Aperture Radar 
(SAR) images as long, slightly curved, internal wave fronts 100-200 km in 
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length and separated by 50-80 km. When these isolated internal wave fronts 
sweep across Dongsha Atoll in the middle of the South China Sea, they 
are perturbed and subsequently re-radiated as two sets of wave front trains 
propagating westward, each with greater curvature than the incoming wave 
front. The SAR images of this process are shown in Figure [1} The transverse 
interactions of the re-radiated wave fronts are so strong that they reconnect, 
rather than pass through each other. This reconnection is shown in Figure 
|2j The reconnection phenomenon for nonlinear internal wave fronts in the 
ocean indicates transfer of momentum and is one of the primary motivations 
of the present study. Internal wave interactions have been well studied in one 
dimension, often by using the weakly nonlinear Boussinesq approximation, 
usually resulting in a variant of the Korteweg-de Vries (KdV) equation and its 
soliton solutions [5T] . However, the complex wave front interactions shown in 
Figures [T] and [2] are plainly two dimensional. Moreover, their reconnection is 
not captured by the 2d extension of KdV for weakly nonlinear waves with slow 
transverse variations, known as the Kadomsev-Petviashvili (KP) equation. 
(The KP equation assumes weak gradients in the direction transverse to 
propagation. However, this assumption does not hold during the wavefront 
reconnections we seek to model.) Thus, we begin our study of internal wave 
front reconnection by developing a new set of equations that extends, to 
multilayer fluids in two dimensions, the Su-Gardner or Green-Naghdi [TH] 
equations for fully nonlinear waves on the free surface of a single-layer fluid 
in one dimension. Related earlier derivations of one dimensional equations 
for multilayer fluids appear in [HI |3Hl EH] . 

Our approach to developing these equations takes advantage of the Euler- 
Poincare (EP) theory of Hamilton's variational principle for ideal fluids, in 
the Eulerian representation [23] , expressed as 55* = with S = J i{u) dt for 
Eulerian horizontal fluid velocity u. The advantage of the EP approach for 
our present purpose is that it provides a hierarchy of equations at various 
levels of approximation that preserves the EP mathematical structure of the 
highest level theory. This approach enables us to strike the appropriate 
balance between accuracy and computational tractability, by choosing the 
appropriate level in the hierarchy of approximations. The symmetries of the 
EP variational principle at each level then endow its resulting evolutionary 
equations with conservation laws via the Kelvin-Noether theorem [23]. For 
example, Kelvin's circulation theorem for the multilayer fluid theory in 2d 
leads to local conservation of potential vorticity (PV) in each fluid layer, even 
though these layers are strongly coupled. 



D. D. Holm k M. F. Staley 



Singular Wave Fronts 



9 





Figure 1: Synthetic Aperture Radar (SAR) images of the South China Sea 
surface taken from the Space Shuttle show a sequence of large amplitude 
internal waves at basin scale. These wave fronts are moving westward after 
being created by tides and currents running throTigh the Luzon Strait on the 
eastern side of the basin. These wave fronts encounter the Dongsha Atoll in 
the center of the basin and re-emerge after colliding with it. 



D. D. Holm k M. F. Staley 



Singular Wave Fronts 



10 




Figure 2: The transverse interactions of the wave fronts re-emerging from 
their encounter with the Dongsha Atoll are so strong that the wave fronts 
reconnect, rather than pass through each other. This indicates transfer of 
momentum. 
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The starting point in the derivation of these equations is the assumption of 
columnar motion; so that the horizontal velocity in each layer is independent 
of the vertical coordinate. We explain how substituting this columnarity into 
the EP Hamilton's principle for multilayer fluids reduces their description 
from 3d to 2d, and results in a hierarchy of equations which arises when a 
sequence of further approximations are introduced into Hamilton's principle 
in the EP framework. 

The common features of the equations in this hierarchy are (1) they ex- 
press Newton's Law for the evolution of momentum in the Eulerian fluid 
representation; (2) their nonlinearity involves both momentum and fluid ve- 
locity; and (3) their momentum and fluid velocity satisfy a nonlocal rela- 
tionship, which must be solved by inversion of an elliptic operator at each 
time step. While keeping these features, we simplify the description by mak- 
ing a series of approximations in Hamilton's principle, until we finally arrive 
at a minimal description which still captures the wave front reconnection 
phenomenon. We then employ this minimal description to investigate and 
classify the wave front interactions analytically and in numerical simulations. 
Our application is for internal wave fronts, which we compare with our nu- 
merical simulations in two dimensions. In a further developmental step, we 
also consider numerical solutions of our description in three dimensions. As 
mentioned earlier, this 3d extension yields singular solutions of EPDiff whose 
momentum is defined on contact surfaces. 

The EPDiff equation and its weak wavefront solutions. The EP 

equation at which we eventually arrive by this sequence of variational approx- 
imations is the following evolutionary integral-partial differential equation, 
expressed in vector form as [28l 122] , 



for which we assume periodic boundary conditions. This is the EPDiff equa- 
tion, for "Euler-Poincare equation on the diffeomorphisms," where m and u 
are vectors, ^^m denotes the time derivative of m, and a is a constant param- 
eter with dimensions of length. EPDiff is the n— dimensional version of the 
Id shallow water equation introduced in [6]. EPDiff arises from a variational 
principle based on the fundamental dynamical properties of fluids. Math- 
ematically, EPDiff describes geodesic motion on the diffeomorphism group 




with m = u — a^Au, (1) 
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with respect to ||u||^i, the norm of the fluid velocity, which is the kinetic 
energy of this vertically averaged model fluid [23]. 

The EPDiff equation Q may be written equivalently as differential Rie- 
mann invariance; namely, as invariance of the momentum one-form density 
along the fluid velocity characteristics. That is, 

-f- (m ■ (ix (g) dVol) = , along ^ = u = G * m , (2) 
at dt 

where u = G * m denotes convolution of the momentum m with the Green's 
function G to produce the fluid velocity u. This convolution is the elliptic- 
solve step in the determination of the fluid's velocity from its momentum for 
this class of equations. This property is the nonlocal relationship mentioned 
above. In the particular case of EPDiff, we use the elliptic Helmholtz operator 
in m = u — a^Au, with length scale a. The Helmholtz operator relationship 
between velocity u and momentum m is derived in Section [5] for strongly 
nonlinear columnar motion of shallow water. 

The relation of EPDiff to Id soliton equations. In Id, EPDiff becomes 

dtm + mdxU + dxijnu) = , with m = u — a^d^u . 

This is the dispersionless limit of the Camassa-Holm (CH) equation 

dtm + mdxU + dx{mu) = —Cod^u — Td^u , with m = u — a^dlu , 

which is a completely integrable soliton equation |i6j. The CH equation re- 
duces to the familiar KdV equation when — )■ 0. 
EPDiff in Id with viscosity v is expressed as 

dtm + mdxU + dx{mu) = ud^m , with m = u — a^d^u . 

When — 7- this becomes the well known Burgers equation, 

3 

dtu + dxi-u^) = vdlu . 

The Burgers equation is hyperbolic and supports weak solutions (shocks) in 
the limit — ?■ 0. 

EPDiff is a new departure for nonlinear hyperbolic equations. It is con- 
servative in the absence of viscosity and it may then be written in Riemann 
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invariant form However, it is nonlocal, because obtaining its character- 
istic velocity u from its momentum m by inverting the Helmholtz operator 
via u = G * m requires an elliptic solve at each time step when a ^ 0. Our 
analysis |28l |22] shows that EPDiff has interesting and unusual solution prop- 
erties which are only beginning to be studied. The properties of its solutions 
will provide great challenges for analysis and numerics. 



Singular momentum solutions of EPDiff. For example, EPDiff has 
weak singular momentum solutions that are expressed as [281 122] , 

TV 

m(x,t) = 5^/ Pait,Sa)6{^-Qait,Sa))dSa, (3) 
a=l -^^a 

where Sa is a Lagrangian coordinate defined along a set of curves in the 
plane by the equations x = Qa(t,Sa) supported on the delta functions in 
the EPDiff solution Q. Thus, the singular momentum solutions of EPDiff 
are vector valued curves representing evolving wave fronts defined by the 
Lagrange-to-Euler map ^ for their momentum. 

The Green's function for the Helmholtz operator relates the fluid velocity 
to the momentum in EPDiff ([T]). Thus, substituting the defining relation 
u = G * m into the singular momentum solution ([s]) yields the velocity 
representation for the wave fronts as another superposition of integrals, 

N 

u(x,t) = ^/ P,{t,S,)G{^,Qa{t,S,))dSa. (4) 

a=l '^•S'o 

The Green's function G for the second order Helmholtz operator in this 
expression has a discontinuity in slope across each Lagrangian curve moving 
with the velocity of the flow. Being discontinuities in the gradient of velocity 
that move along with the flow, these singular solutions for the velocity at the 
wave fronts are classified as contact discontinuities in the fluid. 

Thus, the weak solutions of EPDiff are contacts, rather than shocks. 
This behavior challenges numerics and leads to a wide range of potential 
physical applications. In addition to describing internal wave fronts, the 
EPDiff equation describes the interaction of contacts in a variety of other 
situations ranging from solitons |6] to turbulence [9l |T6] to computational 
anatomy [271 SB] . The EPDiff fluid dynamics equation in ([T]) describes the 
limiting (pressureless) cases of all of these applications. 
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Relation between contact solutions of EPDiff and solitons. The 

EPDiff singular solutions Q for the velocities of the internal wave fronts 
represent the third of the three known types of fluid singularities: shocks, 
vortices and contacts. The key feature of these contacts is that they carry mo- 
mentum; so the wave front interactions they represent are collisions, in which 
momentum is exchanged. This is very reminiscent of the soliton paradigm in 
one dimension. Indeed, in one dimension the sing ular solutions Q of EPDiff 
are true solitons that undergo elastic collisions and are solvable by the inverse 
scattering transform for an isospectral eigenvalue problem [6j. 

Applications of EPDiff in turbulence. Lagrangian averaging is a 
promising approach in turbulence closure modeling. This approach has 
the advantage that Lagrangian averaging commutes with the advective time 
derivative in fluid mechanics. Thus, Lagrangian averaging preserves the vor- 
ticity stretching process in the resulting approximate equations. For exam- 
ple, after Lagrangian averaging the Navier-Stokes equation and using Tay- 
lor's hypothesis that the fluctuations move with the mean flow, one flnds the 
following turbulence closure model [U] , 

-^v + u ■ Vv + Vu^ . V + V» = z/Av + F , (5) 
ot 

with V = u — a^Au and V • u = . 

These are the equations of the Lagrangian averaged Navier-Stokes alpha 
(LANS-alpha) model of turbulence, whose properties are reviewed, for exam- 
ple, in [16]. EPDiff ([T]) is recovered from the LANS-alpha equations ([s]) when 
the constraint of incompressibility (V ■ u = 0) is relaxed (so that pressure 
gradient Vp may be dropped), and the viscous and forcing terms on the right 
hand side are absent. Thus, these equations possess the following analogy: 
EPDiff is to LANS-alpha, as Burgers is to Navier-Stokes. Namely, Burgers 
is a simplifled model of compressible Navier-Stokes which allows shocks as 
singular solutions of its initial value problem, when F is absent and — )• 0; 
and EPDiff is a simplifled model of compressible LANS-alpha which allows 
contact discontinuities as singular solutions of its initial value problem, when 
F is absent and z/ — )■ 0. The reduction of the singular solutions from Burgers 
shocks to EPDiff contacts is a result of the Lagrangian averaging process, 
which tempers the nonlinearity in the Navier-Stokes equations. 
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Applications of EPDiff in computational anatomy. Applications of 
EPDiff ([T]) in computational image science (CIS) focus on computational 
anatomy, based on pattern matching algorithms and smooth morphing of 
planar figure outlines, called templates. For a review of this approach, see 
[lU] . An interesting objective of this CIS application is to reconstruct a 3d 
map of the brain in the spatial region between two parallel 2d PET Scans. In 
this application, the figure outlines are contact curves, which form a finite- 
dimensional invariant manifold of EPDiff, as in equation ([s]). A 3d brain map 
is constructed by treating two 2d PET scans in parallel planes as initial and 
final conditions, and then fiowing by EPDiff to interpolate in the 3d region 
between them as an optimization problem. The 2d solution for the contact 
curves evolving by EPDiff smoothly reconstructs the outlines of the PET 
scan images on any parallel plane between them. The EPDiff contact inter- 
actions also allow reconnection of planar outlines, corresponding to changes 
of topology in planar sections of the 3d object being imaged [27] . The rep- 
resentation ^ of the singular momentum solutions of EPDiff encodes the 
image contours for this application. 

Importantly, the momentum representation of the image contours ^ is 
complete and nonredundant (one-to-one). Another advantage of the momen- 
tum representation of image contours as singular solutions of EPDiff is that 
this representation is linear in nature, being dual to the velocity vectors. 
Thus, linear combinations of either velocity fields, or momenta are mean- 
ingful mathematically and physically, provided they are applied to the same 
template [27]. This means the average of a collection of momenta, of their 
principal components, or time derivatives of momenta at a fixed template are 
all well-defined quantities. The linearity the momentum representation also 
allows for example the statistics of an ensemble of images to be analyzed, 
or the results of adding noise to the image outlines to be computed using 
EPDiff. All of these advantages of the singular momentum representation 
for 2d images also apply in the corresponding representation of 3d images. 
Evolution of singular momentum solutions of EPDiff as surfaces in 3d cor- 
responding in medical imaging to growth, or changes in 3d shape with time. 
The last part of this paper will deal with the computation of 3d solutions of 
EPDiff. 

Summary of the analytical approach. The rich array of possible ap- 
plications of EPDiff motivates our investigation of its initial value problems 
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in Id, 2d and 3d. The paradigm raised in our analysis of internal wave front 
interactions is the generalization of soliton collisions from Id to 2d. The 
class of shallow water equations we derive and study in 2d has a sequence of 
simplifications that yield two dimensional generalizations of all the weakly 
nonlinear shallow water equations that have historically been used to study 
solitons and solitary waves in one dimension. We first argue that by using 
the highest, most accurate level of the fully nonlinear description in 2d, one 
should be able to reproduce the observed phenomenon of wave front recon- 
nection observed for internal waves in the South China Sea. We then simplify 
the description by stages until we finally arrive at EPDiff Q in 2d, while 
preserving what we believe is the key feature responsible for this wave front 
reconnection phenomenon - its momentum transfer during collisions. Be- 
cause of its other potential applications, particularly its linear encoding of 
information for imaging science as momentum, we also consider the initial 
value problem for 3d singular solutions of EPDiff. Its imaging science appli- 
cations are envisioned as optimization problems. However, the first step in 
understanding the role of momentum exchange for EPDiff in imaging science 
is the solution of its initial value problem for the emergence of its singular 
momentum solutions in 2d and 3d. 

Brief sketch of the paper. The first few sections of this paper motivate 
using the much simpler equation EPDiff ([T]), whose properties we study ana- 
lytically and numerically in later sections, as a simplified but realistic model 
of wave front interactions of nonlinear internal waves in the ocean. The last 
few sections present numerical results for EPDiff in both 2d and 3d, which 
we now preview. 

3 Preview of numerical results 

3.1 CFD methods for fluid singularities 

Shocks, vortices and contacts are a compressible fluid's singular nonlinear 
responses to strong applied forces. Capturing these singular responses accu- 
rately has always been the grand challenge of computational fluid dynam- 
ics (CFD). Shocks move through the flow and carry inertia. Vortices move 
with the flow, but they have no inertia. Contacts are discontinuities in the 
derivatives of velocity and density that, like vortices, move with the fluid 
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flow, and which do have inertia. Consequently, it makes dynamical sense to 
speak of momentum exchange in contact-contact collisions. Today, various 
CFD methods exist that accurately capture shocks and vortices. However, 
considerably less is known about designing numerical methods for capturing 
contacts and characterizing their nonlinear interactions, especially in higher 
dimensions, when vorticity is present. 

One might consider Lagrangian numerical methods as an obvious ap- 
proach, because contact discontinuities move with the flow. However, for 
head-on contact collisions, parallel studies showed that the Lagrangian ap- 
proach suffers in comparison to compatible differencing algorithms (CDAs) 
that we apply here in the Eulerian framework. In particular, the head-on 
collision of oppositely moving contacts produces an elastic bounce seen in 
the weak solution of EPDiff as a mutual annihilation and recreation of gen- 
eralized functions, but only the annihilation was captured well by Lagrangian 
methods in the parallel studies [37j. In the future, we may also consider de- 
veloping new methods for EPDiff based on discrete exterior calculus (DEC) 

[iniEH]. 

3.2 Numerical approach 

The second half of this paper numerically investigates and classifies the 
EPDiff dynamics of contact interactions in various cases of its initial value 
problem. This is accomplished by applying compatible differencing algo- 
rithms to EPDiff in both 2d and 3d. One reason for choosing CDAs as the 
preferred simulation method is that the EPDiff equation ([T]) may be rewritten 
as 



This equation contains the divergence, gradient and curl operators, whose 
identities must be treated properly for accurate computations. Preserving 
these vector calculus identities is the basis of CDAs and is the first funda- 
mental property of discrete exterior calculus. Namely, these identities form 
the discrete analog of the identity (P = for the exterior derivative. As men- 
tioned earlier, EPDiff represents geodesic motion, which naturally involves 
exterior calculus and variational principles. A framework for designing DEC 
methods with promising potential for simulations of EPDiff has been devel- 
oped and advanced recently in [191 135] . 
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3.3 Interaction dynamics of contacts 

This paper takes advantage of recent developments in compatible differencing 
algorithms by applying CDAs to new classes of problems involving contact- 
contact interaction phenomena. For the report of a recent conference on 
CDAs, see [52] • We characterize the emergence of contacts from smooth 
initial velocity distributions, and describe their subsequent evolution, prop- 
agation and interaction dynamics. We address dynamical issues for basic 
interactions among contacts in one, two and three dimensions, as follows. 

Id. An integrable shallow water equation whose peaked soliton solutions 
emerge from a smooth spatially confined initial conditions for velocity. 

2d. Interacting contact curve segments in the plane: trains of contact 
curves emerging from an initially continuous fluid velocity distribution, prop- 
agating, and interacting nonlinearly through fundamental collision rules. The 
2d collision rules for singular solutions of EPDiff are elucidated by plotting 
Id linear sections through the 2d solutions which show their spatial profiles 
in various directions. The solution behavior along these sections shows the 
same elastic collision properties and momentum exchange as seen for the Id 
peaked soliton solutions (peakons) in Id for the dispersionless CH equation. 

3d. Interacting contact surfaces in space: sheaves of contact surfaces 
are shown emerging from an initially continuous fluid velocity distribution. 
We investigate their propagation and interaction by plotting level surfaces of 
speed, as well as plotting 2d planar slices through these surfaces. 

More specifically, we examine dynamical issues for the following basic 
interactions among the contacts, which are weak solutions of EPDiff, in one, 
two, and three dimensions. 

Id. In Id, contacts move as points on a line, and EPDiff reduces to the 
well known dispersionless Camassa-Holm (dCH) equation for shallow water 
waves, whose contacts are solitons with sharp peaks, or peakons, as shown 
in the time evolution in Figure |3j The A^— peakon problem is known to be 
completely integrable [6]. 

2d. In 2d, contacts move as segments of curves in the plane. These EPDiff 
contacts correspond to oceanic internal waves in one application, or to pla- 
nar image outlines in another. We shall address the following specific 2d 
scenarios. 
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Evolution and interactions of EPDiff contact segments which are 
initially straight line segments. Contacts move at the local fluid ve- 
locity; so a contact segment must terminate with zero velocity. Hence, an 
initially straight-line contact segment does not remain straight. Instead, it 
evolves under EPDiff into a contact curve segment whose length increases, 
as shown in Figure [ij The speed, |u|, is displayed as colors in this figure, 
while the bottom panels shows profiles in the east (black), north (dotted 
red), northeast (green), and southeast (dotted blue) directions through the 
center of the box. An arrow shows the initial direction of u. See Section [9] 
for a fuller description of the contents of the figures. 

The transverse collision of two contact curve segments may result in re- 
connection (merger, or melding) of the curves, which is admitted by the rapid 
evolution of EPDiff along the directions tangential to the curves. In Figure 
|5} note the striking similarity between the observed behavior, particularly in 
the second frame, and the behavior seen in the earlier image of ocean internal 
waves in Figure |2} We shall investigate the "collision rules" under which such 
reconnections of wave fronts occur. 

We will also consider offset collisions of initially parallel contact segments, 
as shown in Figure |6} In that figure, notice the recreation (after annihilation) 
of the portions of the contact curves in the regions where a head-on collision 
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Figure 4: An initially straight, rightward moving peakon segment balloons 
outward and stretches as it expands in 2d, but keeps its integrity as a single 
wave front. 



n 


5 1 















Figure 5: A skew overtaking collision of two peakon segments shows reconnec- 
tion in which the wavefronts merge repeatedly upon intersecting transversely. 
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Figure 6: Head-on collisions of peakon segments first annihilate and then 
recreate their wavefronts. 



takes place. 

Emergence of contact curve segments from representative classes 
of smooth initial fluid velocity distributions and their subsequent 
interactions. For example, consider an initial velocity distribution whose 
magnitude (speed) is a circular Gaussian ring (annulus) and whose direction 
is chosen in the following ways. 

Uniform translation (broken angular symmetry). This simulation, shown 
in Figure [7} demonstrates that the initial value problem for EPDiff tends to 
produce only contact solutions. It also shows the differences in their propa- 
gation under convergent (left half of annulus) and divergent (right half of an- 
nulus) geometry, and illustrates the exchange of momentum in an overtaking 
collision. The collision as the rightmost contact curve is overtaken from the 
left transfers some of the overtaking curve's momentum forward and causes 
the rightmost curve to bulge slightly, as seen in the last two frames. (Thus, 
while momentum is transferred in this collision, the overtaking contact curve 
does not "pass through" the rightmost contact curve.) 

Uniform rotation (preserving angular symmetry). Figure |8] illustrates how 
the radial and azimuthal components of EPDiff are coupled for radially sym- 
metric solutions. The figure also shows collapse to the origin and reflection 
back outward, which presents an extreme test of the accuracy of our Carte- 
sian numerical algorithm. The radial collapse and bounce can be computed 
semi-analytically for comparison, |26j . 
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Figure 7: A rightward moving circular Gaussian ring of initial velocity breaks 
into divergent and convergent contact wave fronts. 




Figure 8: An initially rotating circular Gaussian ring of velocity couples its 
angular motion to the radial motion of contact wavefronts which propagate 
both inward and outward. The collapsing circular wavefronts reflect from 
the center of symmetry. 
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Figure 9: An initially circular Gaussian ring of velocity is diverging along one 
diagonal and converging along the other. This initial condition breaks into 
interacting contact wavefronts which propagate outward along one diagonal 
and converge inward along the other to annihilate and re-emerge, leaving 
behind a complex mixing flow in the center. 



M-fold discrete angular symmetry. Formation of contacts in the outward 
divergent part of the flow, and the emergence, annihilation, and recreation 
of contacts collapsing to the diagonal, are all illustrated in Figure [9j This 
process is typically followed by very complex patterns of mixing flow involving 
peakon profiles along each Id section. 



3d. In 3d, the initial value problem for EPDiff produces contact surfaces 
moving through space. The interactions (collisions) of these contact surfaces 
produce very complex but coherent patterns. In the figures we discuss here, 
the left, back, and bottom panels of each box show 2d planar slices through 
the center x (southeast), y (northeast), and z (north) values of |u|. The 
corresponding planar slices through the level surfaces of |u| are shown on 
the side panels of each box. 2d slices that form a plane of symmetry in 3d 
invariably reflect the behavior of the corresponding 2d problem. We will 
address the following scenarios in 3d. 



Evolution and interactions of contact surfaces that are initially disc- 
shaped velocity distributions. Initially fiat contact surfaces evolve by 



ballooning into curved contact surfaces, as shown in Figure 10 



Because EPDiff is isotropic, it allows equally rapid evolutions in the di- 
rections normal and tangential to the contact surfaces. In particular, its 
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Figure 10: In 3d, an initially rightward moving disc with exponential trans- 
verse velocity profile (of correct width alpha) balloons outward as a contact 
wavefront which retains its integrity. 




Figure 11: The skew collision of two initially disc-shaped contact wavefronts 
shows reconnection in 3d. 



tangential evolution admits significant stretching. It also allows reconnec- 
tion of contact surfaces which intersect transversely as they collide. The 
reconnection is caused by transfer of momentum. This phenomenon is illus- 



trated in Figure [TT| We shall investigate the "collision rules" under which 
such reconnections of surfaces occur. 

An offset collision of initially parallel contact surfaces is shown in Figure 



12 As in 2d, we observe the recreation (after annihilation) of a contact 
surface after a head-on collision takes place. Notice also how the outer edges 
(away from the head-on collision) of the contact surfaces merge to form a 
ring. 



Emergence of contact curve surfaces from representative classes 
of smooth initial fluid velocity distributions and their subsequent 
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Figure 12: The head-on colhsion of two initially disc-shaped contact wave- 
fronts shows annihilation and reconnection in 3d. 




Figure 13: A rightward moving spherical Gaussian shell of initial velocity 
breaks into divergent and convergent contact wave fronts. 



interactions. For example, consider an initial velocity distribution whose 
magnitude is a radially symmetric Gaussian shell (spherical annulus) and 
whose direction is chosen in the following ways. 

Uniform translation (broken angular symmetry). The simulation shown 
in Figure 13, as in the corresponding figure for 2d, demonstrates that EPDiff 
tends to produce only contact solutions, and illustrates the differences in their 
propagation under convergent (left half of spherical annulus) and divergent 
(right half of spherical annulus) geometry. 

Uniform rotation (preserving angular symmetry about the vertical direc- 
tion). The frames in Figure [M] show collapse of the spherical contact surfaces 
to the origin, and their reflection back outward. This shows the coupling be- 
tween angular and radial motion. 

Two-fold discrete angular symmetry. Figure 15 illustrates the formation 
of contacts for the outward divergent part of the flow and the emergence, 
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Figure 14: An initially spherical Gaussian shell of velocity rotating about the 
vertical axis expands outward and collapses inward with cylindrical symme- 
try, as it breaks into divergent and convergent contact wave fronts. 




Figure 15: An initially spherical Gaussian shell of velocity diverging in one 
diagonal vertical plane and converging in the other breaks into interacting 
contact surfaces. 



annihilation and recreation of contacts collapsing along lines of discrete sym- 
metry. The 2d slice through the horizontal midplane (which is projected on 
the bottom panel) shows 2d behavior similar to that in Figure [9] 

4 History of modeling internal wave fronts 

Modeling observed internal waves propagating over real topography requires 
a fully two-dimensional, multilayer description. Strong, complex, two dimen- 
sional wave-wave and wave-topography interactions have long been observed, 
for example, in internal waves propagating through the Strait of Gibraltar 
|53j . Even more complex interactions were recently observed from the Space 
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Shuttle for internal wave trains propagating in the vicinity of Dongsha Island 
in the South China Sea by Liu et al. [lO]; see Figures [l] and [2] These interac- 
tions produce remarkable nonlinear phenomena. In particular, they produce 
wave front reconnection, as well as the diffraction and refraction expected of 
large amplitude internal waves interacting with bathymetry and boundaries 
such as straits, coasts, shoals, islands and atolls. 

Equations for strongly nonlinear dispersive waves on the free surface of 
a single layer of incompressible fluid with topography are usually attributed 
to Green and Naghdi [ISj, although the same equations were derived earlier 
by Su and Gardner [HU]. These equations for single layer columnar motion 
(ILCM) generalize the Boussinesq family of shallow- water equations to al- 
low for strong nonlinearity. For finite wave amplitudes, the ILCM equations 
capture strongly nonlinear effects, such as the upstream emission of solitary 
waves and the downstream surface depression and oscillations due to flow 
over a obstacle at Froude numbers greater than unity |l3l [39]. References 

[55] show that numerical simulations of the ILCM equations tend to be 
faithful representations of the corresponding Euler solutions, provides that 
wave breaking does not occur. The works of Choi and Camassa [lOl [H] ex- 
tended the ILCM description to a two-layer fluid with fixed horizontal upper 
and lower boundaries, including the cases of two thin layers and of a thin 
layer over an infinitely deep layer. The CC equations admit bi-directional so- 
lutions and they may be derived from a variational principle by following the 
averaged Lagrangian methods pioneered in Whitham [51]. (See also Miles 
and Salmon who derived the ILCM equations using Hamilton's prin- 
ciple in the material representation.) As a consequence of their variational 
principle, the CC equations possess conservation laws for mass, momentum, 
and energy. The same equations were derived earlier and analyzed for well- 
posedness in Liska, Margolin and Wendroff [381 EH] by substituting directly 
in the governing equations the columnar motion assumption that the ver- 
tical velocity is a linear function of vertical coordinate. Choi and Camassa 

showed explicitly, by assuming first weak nonlinearity and then unidirec- 
tional wave propagation, that these equations recover all the known weakly 
nonlinear evolution equations for interfacial waves. The well-posedness is- 
sue identified by Liska et al. [38l [39] for these equations is similar to the 
steepening lemma result proved for the CH equation in [6]. Namely, an ini- 
tial velocity distribution possessing an inflection point of negative slope will 
develop a vertical slope in finite time. This loss of well-posedness in finite 
time is part of the mechanism for the formation of the singular solutions we 
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seek to investigate here. These singular solutions dominate the initial value 
problems we study and one should expect to deviate from standard concepts 
of well-posedness in studying their emergence from smooth initial conditions. 

In the next few sections of this paper, we shall use the Euler-Poincare 
(EP) variational principle for continuum motion in the spatial, or Eulerian, 
representation [23] to extend the CC equations by deriving multilayer colum- 
nar motion (MLCM) equations. MLCM describes strongly nonlinear inter- 
nal waves propagating on the interfaces of layer-stratified incompressible fluid 
driven by gravity, moving over topography and possessing a free surface. The 
MLCM equations are derived here by imposing columnar motion in the EP 
variational principle for the Euler equations of a multilayer incompressible 
fluid with a free surface and variable topography. In one dimension, when 
the free surface and bathymetry are neglected, the MLCM equations recover 
the CC equations [11]. 

Strictly speaking, the CC equations are for vertically averaged horizon- 
tal velocities, not simply for columnar motion. Thus, they need not have 
emerged from an action principle for multilayer columnar motion. The ver- 
tically averaged horizontal velocities do undergo columnar motion, but this 
alone is insufficient to expect a priori that the CC equations would arise 
by substituting columnarity into the action principle. However, Choi and 
Camassa noticed a posteriori that their equations do satisfy a variational 
principle [TT] . We have used their observation, combined with the suggestive 
results of Liska, Margolin and Wendroff [3HI EH] and of Miles and Salmon [15] 
to extend the CC equations to the more general MLCM situation required for 
modeling the large scale internal wave interactions observed in the Gibraltar 
Strait and the South China Sea. 

The new MLCM equations (20) derived here generalize the CC equations 
for the interfacial motion between two layers with fixed horizontal upper and 
lower boundaries, to allow for waves at the top free surface as observed by 
the Space Shuttle, to include an arbitrary number of fiuid layers, and to ac- 
count for variable bottom topography. These equations provide fundamental 
insight into how topography and the multiple layers of a stably stratified 
incompressible fiuid are coupled by nonlinear, nonhydrostatic, dynamical ef- 
fects. This dynamical coupling arises in addition to the more familiar mul- 
tilayer hydrostatic effects. We derive these nonhydrostatic equations and 
their natural boundary conditions from the Euler-Poincare variational prin- 
ciple for continuum motion of Holm, Marsden and Ratiu [23] by including 
in the fluid Lagrangian the kinetic energy due to vertical oscillations in the 
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columnar motion approximation. Their traveling waves, stability properties, 
well-posedness, weakly nonlinear aspects, relations to other approximations 
and numerical simulations will be addressed elsewhere. 

After deriving the highest level nonhydrostatic MLCM equations for mul- 
tilayer internal waves, we make a series of Boussinesq-like approximations for 
weakly nonlinear waves that allows comparison with previous work and which 
eventually results in a minimal description, EPDiff ([T]). We then discuss some 
of the geometrical properties and singular solutions of EPDiff. Finally, we 
describe our numerical approach and present a large suite of 2d and 3d sim- 
ulation results for EPDiff. 



5 Nonhydrostatic multilayer columnar mo- 
tion (MLCM) equations 



Consider a multilayer fluid consisting of immiscible layers moving under 
the constant vertical acceleration of gravity. Regard the top layer (whose 
rest position is z = 0) as first and the bottom layer as last (iVth), so that i 
increases with the depth of the layer. The i—th layer has constant density pi, 
horizontal velocity Ui{x,y, z,t), vertical velocity Wi{x,y, z,t), upper surface 
at z = hi{x, y, t) and lower surface at z = y, t), so the layer thicknesses 

are = hi — /ij+i, with i = 1, . . . ,N. The A^th layer (on the bottom) has 
density p^, horizontal velocity ujy, upper surface at z = h]\r{x,y,t) and 
lower surface at z = /iat+i = —b{x,y), the fixed bottom topography. We 
shall assume the multilayer fluid is stably stratified, so that pi < Pi+i, for 
i = 1, . . . , N (density increases downward). 

5.1 Implications of the columnar motion ansatz 

The Lagrangian in Hamilton's principle for a multilayer fluid is the difference 
of its kinetic and potential energies. 



i = dxdy 




i=l 



(7) 



One may also include the effects of Coriolis force due to rotation by adding 
the term J2iLi Dill{x, y) ■Ui + DiS{x, y)wi to the Lagrangian density, where 
2VL = curl (R(x, y) + S{x, y) z) is twice the rotation vector. 
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We note that incompressibility V ■ Uj = —dzWi relates the horizontal 
and vertical velocities in each layer. When one also assumes the columnar 
motion ansatz, namely, 

then incompressibility also implies linear dependence of the vertical velocity 
on the vertical coordinate in each layer, as 

Wi = -zV ■ Uj + dthi+i + V ■ (/ij+iUj) , (9) 

for hi^i < z < hi. The thickness Di = hi — hi^i of the z— th layer obeys a 
continuity equation, 

dtDi + V ■ (Au.) = . (10) 



Hence, the volume of each constant density layer is conserved and the layer 

D 



thicknesses remain positive. We also note that /ij+i = — h{x, y) + Yl!j=i+i 



as the sums of differences Dj cancel in pairs and hN+i = — b{x, y). The total 
depth is Y^f^i Dj = hi{x, y, t) + b{x, y). 

As a result of these relations, the vertical velocity in the z— th layer for 
hi+i < z < hi may be expressed in terms of the horizontal velocities, the 
layer thicknesses, and the prescribed bathymetry, b{x,y), as P 

N 

Wi = -zV ■ Ui -V ■ b{x, y)ui- ^ V • Dj (u^ - u^) . (11) 

j=i+i 

We shall substitute this expression obtained from columnarity into the 
Lagrangian ([T]), then perform the vertical integrals and use the Euler- 
Poincare theory to obtain the motion equation in each layer, as the EP 
equation |23] . 

-— + u,-V— + Vuf ■ — + — V-u,-AV— = 0. 12 

Ot OUj OUj OUj OUj oDi 

This procedure will produce the multilayer columnar motion (MLCM) 
equations (20), which are completed by the corresponding continuity equa- 
tion ( 10 ) in each layer. 
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5.2 Variational derivatives of the Lagrangian 



Columnarity and its implied formula (11) for Wi allow the vertical inte- 



grals in the Lagrangian ([t]) to be performed as, 

1 f ^ r 



2 



(13) 



dxdy . 



Perhaps not unexpectedly, the integrated kinetic energy is a layer-thickness- 
weighted sum of squares of the horizontal velocities and their divergences. 
The additional notation is defined as ^4^ = V ■ Uj , so that 

D,Ai = - (d/dt + u, ■ v) A = - D\ = n,+i -Ki , 

and 



TV 



j=i+i 



''^i\z=hi+i 



(14) 



Thus, DiAi is the rate of expansion of the i— th layer, and Bi is the vertical 
velocity at its lower interface. Moreover, the vertical velocity at its upper 
interface is DiA^ + Bi = —fii. The differences of squares in the potential 
energy of (13) may also be written in terms of the layer thicknesses upon 

(15) 



substituting for /ij+i, as 



1 



1 ^ 



j=i+i 



Thus, columnarity (|8]) allows the Lagrangian ([T]) to be expressed solely in 
terms of the horizontal velocities {u}, their (weighted) divergences, and 
the layer thicknesses {-D}. The potential energy of each layer is coupled 
hydrostatically to all the layers beneath it by the last term in (15). We 
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leave the top layer free. For a rigid lid [2] one would add the constraint 
hi = —b{x, y) + J2iLi Di = 0, imposed by a Lagrange multiplier ps (the 
surface pressure) determined by preservation of hi = 0. 

We rearrange the sums by using the following identity, which holds for 
arbitrary quantities Qi and Ri, with i = 1,2, . . . , N, 



N 



N 



N 



j-1 



E h^^ E =E h^E^^ 



(16) 



2 = 1 



j=i+l 



i=l 



Consequently, we find that the variational derivatives of the Lagrangian ([T]) 
under columnarity (fsl) are given by |3] 



N 

dxdy y^pi< 



i=l 



i-1 

l"^ - ghi- — ^pjDjj6Di 



+ 



+ 



i—l i—1 

-{D,A, + B,f + E V ■ C,M, - u, ■ V 5^ C 



5Di 



(17) 



where we have introduced notation for two more linear combinations of the 
interface velocities. 



= l-DiA, + l:Bi = - \{2hi + hi+i] 
o / o 

and Ci = DiGi , with average interface velocity, 

= ^ + Bi = - hhi + hi+i) . 



(18) 



(19) 



5.3 Euler-Poincare motion equation for MLCM 



Substituting the variational derivatives from the formula for 51 in (17) into 
the Euler-Poincare equation (12) yields the desired multilayer columnar mo- 
tion (MLCM) system of equations, which may be manipulated into the simple 
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form, 



dui 



dt 



dt 



dGi 
"dt 

dGi 

~dr 



Vb{x,y) 



(20) 



N 



j=i+i 



The top line of equation (20) recovers the ILCM equation of Su-Gardner [50] 
or Green-Naghdi [18] in each layer. The new terms involve sums on j i 



which couple the layers. The nonhydrostatic contributions on the right side 
are dynamical, with d/dt = {d/dt + Uj ■ V). The hydrostatic pressure in the 
i—th layer is given by pi = pigHi , with Hi as in 6i (17), 



H, = hi + 



1 '"^ 



N 



and hi 



(21) 



The quantities Fi and Gi in (20) are defined in terms of the layer thicknesses 
and velocities by equations ( 18[19 ). The resulting MLCM equations (20) are 
closed by the layer continuity relations (10) for Di. 



The MLCM motion equations (20) reduce to the standard ILCM equa- 
tions [501 HH] in the single layer case, in which there are no sums on j i. 
In one dimension, the MLCM equations generalize equations (3.19-3.22) of 
Choi and Camassa [llj for the interfacial motion between two layers with 
fixed horizontal upper and lower boundaries, to allow for waves at the top 
free surface, to include an arbitrary number of fluid layers and to account for 
variable bottom topography. One recovers the CC equations of [HI [381 EH] 
by specializing to two layers with fixed upper and lower boundaries in one 
dimension. This amounts to setting N = 2 and neglecting terms involving 



Bi, Gi and Gj in equations (17 20). 



Nonhydrostatic contributions to momentum and pressure. In the 



variational formula for 6i (17), we see that the horizontal momentum in the 



columnar motion equation for the z— th layer involves horizontal gradients 
mi = — = '^Cij{{D},b)uj 



(22) 



i-i 
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The last equality defines the symmetric, positive definite operator in rrij = 
X]j=i 6)uj, which depends on the set of layer thicknesses {D} and 

the bottom topography b. In terms of this operator, the multilayer fluid 
Lagrangian ([T]) under columnarity ([s]) becomes 



1 f ^ 

-y 5^ (m, • u, - g{h^ - hl^S)dxdy . (23) 



Because the multilayer Lagrangian (23) depends on the horizontal spatial 



coordinate through the bathymetry b{x,y), the spatial integral of the total 



momentum (22) is not conserved, except in translation-invariant directions 
of b{x,y). The total pressure in the i—th layer, 6i/6Di, also contains non- 
hydrostatic contributions arising from the vertical columnar oscillations and 
velocity differences among the layers, in addition to its hydrostatic compo- 
nent. All of these additional nonhydrostatic contributions to the momentum 
and the pressure gradient arise from the kinetic energy of vertical motion 
and are proportional to the velocities of the interfaces 

Kelvin circulation theorem and potential vorticity. Although the 
layer motions are strongly coupled, each layer has its own Kelvin circulation 
theorem, 

c? /" m,- , , 

dx = 0, (24) 



Jc(ui) A 

where the closed loop c(uj) moves with the horizontal velocity Uj in the i—th 
layer. In addition, each layer also locally conserves its own potential vorticity 
(PV), i.e., 

— f + Uj ■ Vgi = where qi = — z ■ curl — ^ . (25) 

Consequently, one has an infinite set of conservation laws, 

C$ = y Di^{qi)dxdy, for any function $. (26) 

The C$ are Casimirs of the Lie-Poisson Hamiltonian operator for the MLCM 
system and they play a role in classifying the MLCM equilibrium solutions, 
as in EQl ET 
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Both the Kelvin circulation theorem and its associated local conservation 
of potential vorticity follow from the invariance properties of the EP vari- 
ational principle. Namely, the EP variational principle is invariant under 
fluid parcel relabeling that preserves Eulerian quantities. For full details of 
the Kelvin- Noether theorem, see [2S]- The implications of PV conservation 
for multilayer internal waves remain to be investigated. PV conservation is 
a new element of the MLCM equations, which differs fundamentally from 
the standard Charney-Drazin non-acceleration theorem approach for purely 
potential waves [8J, in allowing solutions representing waves, bores and cir- 
culations to coexist and interact nonlinearly. 



Lie-Poisson Hamiltonian structure of MLCM. The Euler-Poincare 
variational formulation implies the Lie-Poisson Hamiltonian structure of 
MLCM, upon Legendre transforming the multilayer Lagrangian i in (23), 
as shown in [23]. Hence, the conserved Hamiltonian for the MLCM equa- 
tions is Hmlcm = / Sill ™i ■ Uidxdy — i. Thus, the MLCM equations 
of motion are expressible in Lie-Poisson Hamiltonian form using the stan- 
dard Lie-Poisson bracket in terms of the momenta {m} and the layer 
depths {-D}, as in [231 EHl 121] ■ As expected, the Casimirs for this Lie-Poisson 
bracket are the potential vorticity functionals C$ = / Di^{qi) dxdy, which 
satisfy {H,C<s,} = for any Hamiltonian H. The corresponding treatment 
of the Lie-Poisson Hamiltonian structure for the single layer ON equations 
is given in j^Hl E] • 



Comparison with alternative 2d averaged shallow water equations. 

In addition to possessing conservation properties for energy, circulation and 
potential vorticity, the EP motion equation (20) takes a simpler form than 
many other depth-integrated motion equations discussed in the literature, 
such as in jHl |T71 112]. Whether MLCM will be as successful in simulating 
internal wave interactions remains to be seen. As in all columnar motion 
equations, a key feature of MLCM is the elliptic operator in (22) relating 
velocity and momentum. Our numerical simulations in section |9 show that a 
weakly nonlinear approximation of the elliptic inversion in the MLCM model 
does capture the characteristic aspects of the internal wave-front collisions 
and reconnections which were observed in the South China Sea by Liu et al. 
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6 Weakly nonlinear limit equations 



6.1 EP derivation of Boussinesq-like equations 

One may rewrite the total Lagrangian ((T3| without approximation as 



1 r ^ 



pi 



i=l 



A|u,f 



dxdy . 



(27) 



The two-dimensional CC equations arise as EP equations from this La- 
grangian upon neglecting its final term, | J piDiBi_iBi dxdy. We in- 
troduce an alternative approximation of the last term, by approximating it 
in the weakly nonlinear limit, as 



1=1 L 

+ |((divAui)' + 3S,_iS, 



dxdy . 



(28) 



where {di : i = 1,2, . . . , N} are a set of constants, differing only slightly 
from the corresponding Di. In this alternative approximation, the La- 
grangian (28) has the following variational derivatives 

di 



lUi = 



_ / di , 1 6Bi 

piDi \ii - — V(div DiUi) + T^r ^ 
V 6 Uj ou 



5\ii 



for the momentum density, where Bi = {di/2)Bi_iBi , and 
1 6i 

~ 2 



Pi SDi 



-\ui\ - gHi - — V(div Aui) + YTT 
6 oJJi 



(29) 



(30) 



for the Bernoulli potential. 

From these variational derivatives, the multilayer EP equations (12) pro- 
duce 



d ^ ^rr did^Di 

-u, + u, ■ Vu, + gVH, + ^V — 

. 1 SB, ^ / 1 5B, 5B, 

= UiX curl—- V TT^i • ^ jjr- 

Di ou,- \Di o\ii oDj 



(31) 
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We call these the "Boussinesq-like" multilayer equations, because the ap- 
proximation DiAi = DjdivUj ~ divDjUj = ~dDi/dt in the Lagrangian (28) 



replaces the strongly nonlinear dispersive term in the MLCM equations with 
Boussinesq-like linear dispersion, as 



-V d 



djDiAi 
dt 



dj,d_ 



V(div DiUi 



di^d'^Di 

—V 

3 dt^ 



(32) 



This term is linear, by virtue of the continuity equation for each layer thick- 
ness. The completion of this approximation depends on the treatment of the 



Bi terms in equation (31), which are neglected in the CC equations. In the 



approximation that one neglects the Bi terms, the multilayer Boussinesq-like 
equations become 



— u, + u, ■Vui+ gVHi + 

ot 3 ot^ 



0. 



(33) 



These equations are still coupled by the multilayer hydrostatic pressure gra- 



dient, where hydrostatic pressure head is given in equation (21). 

The contributions of the Bi = {di/2)Bi^iBi , terms to the EP equations 



as follows. 



(|29||31|) may be obtained by computing the required variational derivatives, 

di 



N 
N 



-Bi^iBidxdy 



N 



^Pi^{Bi^i 



B.+i)5B, 



i=l 



J2pi^{B.-i + B,+i)Vb-6ui 
i=i 



(34) 



N 
i=l 



i-1 



In the first line, we have used the definition of Bi in equation ( 14 ). In the last 



line, we have dropped boundary terms when integrating by parts and used 
the summation identity (16). These formulas determine the contributions 



of the Bi terms to the EP equations (29 31). In particular, they couple the 



horizontal motion of the i—th layer to the vertical interface velocity of its 
next nearest layer below, since Bi_i + Bi^i = —Ki — Ki+2- Thus, ignoring 
the contributions of the Bi terms corresponds to ignoring interactions among 
next nearest layers in the weakly nonlinear limit. 
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6.2 Other weakly nonlinear EP equations 

We invoke the weakly nonlinear limit to make a further approximation in the 
multilayer Lagrangian (27), by representing its kinetic energy due to vertical 
motion as 



1 1 T.p-^ 



i=l 



dxdy . (35) 



Consequently, its variational derivatives become 



m. 



5u,- 



p,A(ui-^V(AdivUi; 



for the momentum density, and 
I 51 1, 



Pi SDi 



u,- 



6 



(divui)^ 



for the Bernoulli potential. The corresponding EP equations are, 
d 



dt 



Uj + Uj ■ Vui + gVHi 

d d^ d^ 
- ^vk- V(Adiv u,) + Ui X curl V(Adiv u^ 

at 6D.; 6Di 



3 



(divu,)2 + -^u, ■ V(Adivu,)) . 



In the weakly nonlinear limit, one neglects terms in V-Dj to find 

d_ 
dt 



^Mi - ^ V(div Ui) ) + Ui ■ Vui + gVHi 



- y (divu,)' + u, ■ V(divu,)) + 0(VA) . 



(36) 



(37) 



(38) 



(39) 



Upon neglecting terms of order O(V-Di), choosing the center of volume frame 
of reference with -DjUi = for N = 2, and specializing to one dimension, 
one recovers the weakly nonlinear limit of the CC equations in [H]. These 
weakly nonlinear limit equations were shown in [llj to recover all of the 
various Boussinesq approximations for one dimensional shallow water theory. 
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Remark on inversion of the second order elliptic operator. As with 
all the equations in the GN family, the evolution equation (39) requires inver- 
sion of a second order operator in solving for fluid velocity Uj from momentum 
density nij at each time step. A simpliflcation occurs in the weakly nonlin- 



ear approximation, because the operator (1 ^-Vdiv) appearing in this 

approximation is independent of the layer thicknesses, Di, which do enter 
at the fully nonlinear level as in equations (22), (29) and (36). The nonlo- 
cality and smoothing associated with inversion of such elliptic operators is a 
characteristic feature of the entire family of GN equations. 

The remainder of this paper is devoted to using the EP theory in char- 
acterizing the two-dimensional effects of this elliptic operator inversion on 
the solutions of equations in the GN family of shallow water equations. To 
focus our attention on this aspect of the investigation, we shall not require 
the effects of potential energy. 



7 Kinetic energy Lagrangians and the EPDifF 
geodesic equation 

7.1 Kinetic energy Lagrangians 

Neglecting potential energy entirely in the weakly nonlinear limit multilayer 
Lagrangian (35), by setting Di di, gives 

N 

2 /' 



1 f ^ 

2 / Y.P''^' 

i=l 



+ y(divui)' 



dxdy . 



(40) 



The corresponding momentum density involves the second order elliptic op- 
erator, 



m, = 



5ui 



Pidi [vLi 



d; 



V(divui) 



(41) 



Without potential energy, the layers decouple and the EP equation (12) takes 
the same form in every layer, so we may drop the layer index i and write the 
EP equation for the kinetic energy Lagrangian (40) as 



dt 



m + Vu ■ m + m(div u) = , with m = u V(div u) . (42) 
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By design, this equation has no contribution from potential energy. In addi- 
tion, its evolution conserves the kinetic energy. 



|up + ^(divu)^ 



dxdy . 



(43) 



Evolution by kinetic energy in Hamilton's principle results in geodesic mo- 
tion, with respect to the velocity norm provided by the kinetic energy La- 
grangian. 

Reduction to the Camassa-Holm (CH) equation in Id. In one di- 



mension, equation (42) simplifies to 



dm d , , du 
— — h — [mu) + m— 
ot ox ox 



, with m = u — a' 



dx"^ 



(44) 



This is the dispersionless limit of the Camassa-Holm (CH) equation [6], which 
is known to be the equation at quadratic order in the shallow water asymp- 
totic expansion, one full order beyond KdV, whereas KdV appears at linear 
order in this expansion [Til [15]. (The dispersionless limit of the CH 
equation appears because we are ignoring potential energy in this part of our 
investigation.) 



Strengthening the kinetic energy norm to if^. The kinetic energy 
(43) is only part of the H]^ norm of the velocity, defined as 



u 



HI 



|up + (div u)^ + (curl u)' 



dxdy 



lur + a iVul 



dxdy . 



(45) 



Here we assume u is tangential on the boundaries upon integrating by parts 
in Cartesian geometry. We have also simplified the notation slightly by re- 
placing (i^/3 with a^, where a is a length scale. In anticipation that mathe- 
matical analysis will be facilitated by controlling the entire H]^ norm of the 
velocity, we shall choose our kinetic energy Lagrangian to be 



1, 



u 



Hi 



I |2 I 2|VT |2 

u + a Vu 



dxdy . 



(46) 
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The corresponding EP equation (12) is EPDiff, which involves the famihar 
Helmholtz operator in the relation between fluid velocity and momentum 
density, 

d Si 

— m+u-Vm+Vu"^-m+m(divu) = 0, with m = — = u— a^Au. (47) 

ot ou 

We shall assume periodic boundary conditions for the remainder of our inves- 



tigations in this paper. An alternative way of writing EPDiff (47) is 
d 

— m — u X curl m + grad(u ■ m) + m(div u) = , (48) 

which involves the three differential operators curl, gradient and divergence 
in two dimensions. This is the EPDiff equation whose solution behavior for 
the initial value problem is studied in the remainder of the paper. 

7.2 Momentum maps for singular solutions of EPDiff 

Substituting the singular momentum solution formula ([s]) for s G and 
its corresponding velocity ^ into EPDiff, then integrating against a smooth 
test function, implies the following Lagrangian wave front equations, 

dt 

6=1 



Qa{s,t) =J2[ Pb{s',t)G{Qa{s,t),Qb{s',t))ds', (49) 

6=1 ^ 

-P,(s,t) = (Pa(5,t)-P6(3',t))^Q-^G(Q,(s,t),Q,(.',t))rf.'. 



in which summation is explicit on h G 1, 2, ... A^, and there is no sum on 
a. The dot product ■ P;, denotes the inner, or scalar, product of the two 
vectors P^ and P;, in M^. Thus, the momentum solution formula ([s]) yields a 



closed set of integro-partial-differential equations (IPDEs) given by (49) for 
the vector parameters Qa(s, t) and Pa(s, t) with i = 1,2 . . . N . 

Canonical Hamiltonian dynamics of wave fronts in M^. The singular 
momentum solution formula (jsj) is shown to be a momentum map in [22] . 
This fact guarantees the following result. 
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Theorem. The Lagrangian wave front equations (49) are canonical Hamil- 
tonian equations, 



Ft 



Pa{s,t) 



SH 



N 



(50) 



The corresponding Hamiltonian function is, 



N 

HN = \jjY. (Pa(s,t) ■ Pfe(s',t)) G{Cla{s,t)Mb{s\t)) dsds'. 



(51) 



a, 6=1 



Because the solution formula ([3]) is a momentum map the singular solution 
dynamics is collective. That is, the Hamiltonian arises by substituting the 
singular momentum solution formula ^ into the kinetic energy norm ( 46 ) 
and using the delta functions to perform the integrals. Thus, the evolutionary 



IPDE system (49) represents canonically Hamiltonian motion on the space 
of curves in M^. Moreover, this Hamiltonian motion is geodesic with respect 



to the CO- metric given on these curves in (51) by the Green's function G. 
The Hamiltonian Hj^ = ^||P|P in (51) for this motion defines the norm 
||P|| in terms of this co-metric. This momentum map result helps organize 
the theory and provides new avenues of exploration, as suggested in [22] . 
The remainder of this paper, however, will deal with numerical simulations 
which capture the momentum exchange properties of these singular EPDiff 
solutions. 



8 Numerical approach 

Our numerical studies of EPDiff were performed on uniform, logically rectan- 
gular Eulerian grids in 2d and 3d, using the compatible differencing algorithm 



(CDA) described in [31] and sketched in Figure 16 In contrast to our experi- 
ence with Lagrangian methods, our numerics using this CDA have captured 
the elastic bounce expected in head-on collisions with only small distortions 
observed in the recreated contact curves. Future investigations may allow 
us to improve CDAs by developing related variational integrators based on 
additional ideas from discrete exterior calculus (DEC) [T9| [35]. 

In this CDA, scalar and vector quantities are defined at locations that are 
naturally appropriate for the domains and ranges of the discrete divergence, 
gradient, and curl operators. Eight spaces, or grid centerings, include a node 
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Figure 16: This schematic illustrates a particular compatible differenc- 
ing algorithm for quantities defined on uniform, logically rectangular grids. 
Reading left to right, we have nodes, edges, faces, and cells. Nodes and cells 
support scalar- valued functions, while edges and faces support vector- valued 
functions. Divergence, gradient, and curl operators map between nodes, 
edges, faces, and cells. 
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space, at left, which supports scalar functions. An edge space, consisting of 
x-direction edges, ^/-direction edges, and z-direction edges, supports vector 
functions; x components of vectors exists on x-direction edges, etc. A face 
space, consisting of faces perpendicular to the x, y, and z directions, also 
supports vectors. Finally, a cell space supports scalars. 

The set of spaces shown in the figure supports two versions each of diver- 
gence, gradient, and curl, as described in p^. Divergence maps vectors to 
scalars, gradient maps scalars to vectors, and curl maps vectors to vectors. 
Lines between particular spaces in Figure 16 illustrate how the quantities 
defined on different spaces contribute to one another through the various 
operators. 

Node-space scalars map to edge-space vectors through the x, y, and z 
components of the discrete gradient operator, while cell-space scalars map 
to face-space vectors through the x, y, and z components of another discrete 
gradient. The x, y, and z components of edge-space vectors each contribute a 
term {du/dx, du/dy, or du/dz) to a discrete divergence defined on the nodes, 
while the x, y, and z components of face-space vectors each contribute a term 
to the cell-based divergence. Finally, the discrete curl operators map between 
edges and faces in the manner one expects: y and z inputs contribute to x 
outputs, X and z inputs to y outputs, and x and y inputs to z outputs. 

Different components of the same vectors have different numbers of dis- 
crete points in this CDA. If the number of nodes is x M x P, then the 
number of x-direction edges is A^ — 1 x M x P, the number of y-direction 
edges is A^ X M — 1 X P, and the number of z-direction edges isA^xMxP — 1. 
The different numbers of discrete points and their slightly different locations 
must be managed appropriately in any code that uses this scheme. 

At present, our numerics are limited to uniform, logically rectangular Eu- 
lerian grids with periodic boundary conditions. Our results are promising, 
but future work will be needed to continue investigating numerical methods 
and identifying the best candidates for capturing contact behavior on non- 
uniform or unstructured grids and on non-rectangular domains with bound- 
ary conditions other than periodic. For example, when studying the relation 
of EPDiff to internal waves, one might examine the behavior of contact seg- 
ments as they interact with islands or atolls placed into the domain. 

For our 2d and 3d numerical simulations, we advanced the momentum 
m in (48) with an explicit, variable time step Runge-Kutta type predictor- 
corrector. We selected the time step for numerical stability by trial and error, 
while our code selected the time step for numerical accuracy (not to exceed 



D. D. Holm & M. F. Staley 



Singular Wave Fronts 



45 



the time step for numerical stability) according to the following well-known 
formula, 

= 7/^.-1 (52) 



which is used in the following way. At step i of the calculation, we know 
the predicted solution Ui, the corrected solution Ui, and the previous time 
step hi^i. The predictor's order of accuracy is p, while the corrector's order 
of accuracy is p + 1. For our 2d simulations we used p = 4, while we used 
p = 3 for our 3d simulations because this reduced the number of large, 3d 
temporary arrays needed for the calculations, and thereby allowed us to have 
a higher resolution while still using a reasonably accurate time integration 
scheme. 

A new time step is chosen from (52) based on the old time step hi_i 
and the norm of the difference between the current predicted and corrected 
solutions. For both 2d and 3d, we used a relative error tolerance per time 
step of e = 10~^, a safety factor 7 = 0.9, and the L2 norm, || ■ II2. 

Divergence, gradient, and curl were computed at second order accuracy 
according to the CDA outlined above, with the vectors m and u defined 
on the edges (shown as the red, green, and blue spaces in Figure 16). We 
observed that the quality of our numerics did not improve markedly with 
fourth or sixth order operators, and as expected, the higher-order operators 
were somewhat slower to compute. Note that our second order operators, 
in addition to mimicking important properties of their continuum analogs, 
also have greater accuracy than one might expect at second order. This is 
because the staggered nature of the grid allows for a smaller "Ax" in the 
difference computations. For example, du/dx mapping nodes to x-direction 
edges uses the stencil du/dx{i + ^,j,k) = {u{i + l,j,k) — u{i, j, k))/ Ax, 
whereas the analogous computation on a strictly nodal grid is du/ dx{i,j, k) = 
{u{t + k) - u{i - 1, J, k))/{2Ax). 

In two dimensions, one must regard the individual spaces in the schematic 
of Figure 16 as compressed vertically, so that they are fiat. This corresponds 
to having no z component. In this direction edges (red) are identical 

to y-direction faces (yellow), y-direction edges (green) are identical to x- 
direction faces (purple), z-direction edges (dark blue) are identical to nodes 
(brown), and z-direction faces (light blue) are identical to cells (gray). (This 
is not generally true on nonuniform grids.) Even so, the proper treatment of 
quantities defined on the different spaces requires that we regard the different 
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spaces as distinct in 2d. 

In our formulation, the curl operator appearing in (48) is still meaningful 
in 2d. In particular, curl(m) takes the 2d function m, defined on the x- 
direction edges (red) and y-direction edges (green), and maps it to a face- 
space quantity with x and y components of (which are not stored in a 
computer code, of course), and a nonzero z component (light blue in Figure 



16), which is regarded as a vector that is normal to the plane. 

For our 2d simulations, we used a resolution of 1024^ zones, or 1025^ 
nodes. For our 3d simulations, we used 256^ zones (257^ nodes). To invert 
the Helmholtz operator in transforming between m and u, we convolved 
m{x,t) with the Green function in Fourier space. No artificial viscosity or 
other numerical tricks proved to be necessary for our simulations. 



9 Numerical results for EPDiff in 2d 

Upcoming sections describe 2d simulations of evolution under EPDiff for each 
of nine initial velocity distributions. For each simulation we have figures for 
a = a, a = a/2, a = a/A, and a = a/8. Each figure contains six frames, 
showing the initial magnitude (speed) of velocity, |u|, in the upper left frame, 
followed by plots of |u| at future times, reading across and then down. For 
each simulation, the domain is [—1, 1] x [—1,1], x is toward the right, and y 
is toward the top. 

Colors in each frame indicate the magnitude of the velocity, beginning 
with gray for |u| = and ending with white for the maximum value of |u|, as 



shown in the color bar in Figure 17 Maximum values of |u| are determined 
for each frame individually, not over all frames in a figure, so that the colors 
in frames with smaller maximum values of |u| are not washed out. Notice 
how the use of the color black for small |u|, just above gray for |u| = in the 
color scheme, allows us to etch the outlines of the spatially confined velocity 
distributions. 

The transverse profile of the velocity distribution along the horizontal 
midline of each frame is shown as the black (solid) graph in the lower panel 
of each frame, while the red (dotted) graph shows the vertical midline. The 
profiles along the northeast and southeast diagonals are plotted in green 



(solid) and blue (dotted), as sketched in 18 Unlike the colors in the full 
2d plots, which are scaled according to the maximum of |u| in each frame, 
the vertical axis of the profiles in the lower panels are set between and the 
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Figure 17: Color scheme for the magnitude of velocity, |u|, in the upcoming 
figures. 




Figure 18: Locations of Id profiles of |u| shown in the bottom panels in the 
upcoming 2d figures. 



maximum over all frames (in each figure) of the black, red, green, and blue 



profiles. So, for example, in the profiles of Figure 19, |u| as seen in the lower 



panels is seen to decrease as the evolution proceeds, even though the colors 
in the 2d plots always vary between gray (minimum) and white (maximum) 
in order to show maximum detail. 



9.1 Plate 

Each figure we shall discuss compares evolution under EPDiff with various 
values of alpha, starting from the same standard initial velocity distribution. 



which for Figures [T9p2] we call "plate." The velocity in these four 2d plate 
flows is initially rightward. The initial speed in each of these cases is dis- 
tributed along a line segment in 2d, which in 3d will be a disc, or plate. 
Hence, the term plate flow. This initial speed is constant along most of the 
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segment's length, then falls off as a Gaussian at either end. It falls off expo- 
nentially in the transverse directions. The width a of the initial exponential 
profile e~'^'/°" is the same in each figure. We start the EPDiff evolution in 
each figure with this profile, and we vary the parameter a, in the relation 
m = u — a^Au, from problem to problem in fractions of the standard initial 
profile's exponential width a. 

In Figures 19 19, the values of a are, respectively, a, a/2, (t/4, and a/8. 
In Figure 19, the initially straight velocity segment evolves to '"balloon" out- 
ward, while maintaining its transverse peakon profile. In Figures 20 21, the 
initially straight velocity segment in each case evolves into a series of curved 
peakon strips of width a with transverse peakon profiles. The lower panels in 
each case confirm that the velocity profiles in every transverse direction have 
the characteristic exponential peakon shape, e"'"^'/". The first strip to emerge 
travels fastest and each subsequent strip moves slower. Consequently, the dis- 
tances between the strips increases. These peakon strips are curved because 
their endpoints are nearly fixed, while their middle regions are still moving. 
They stretch during their evolution and increase their lengths. Because of 
their peakon cross-sections in velocity profile, each peakon strip corresponds 
to a singular momentum density supported on a curve through the center 
of the strip. The speed varies along each strip according to its height at 
the peak at a given point along the curve. Consequently, the "hotspots" 
appearing in the velocity color scheme as red inhomogeneities along the strip 
are moving faster than the adjoining yellow regions. The extension of these 
hotspots along the peakon strips also indicates order 0(1) stretching. 

One sees kinks near the endpoints of the curved peakon strips that first 
arise at a distance of order a, the matching length in the initial velocity 
profile. Thus, the connection along the segment to the zero speed background 
influences the stretching of the peakon segments over a finite length scale. 
These kinks near the endpoints are more pronounced for the larger values of 
a than for the smaller values, indicating that the length scale a plays a role 
in the stretching process, as well as in the shape of the transverse profile. 
These peakon segments also rotate around their nearly fixed endpoints as 
their expansion and stretching proceeds. 

The black (horizontal) transverse profiles in the lower panels of the fig- 
ures show a steady rightward progression of peakon profiles. The green and 
blue profiles show a similar progression of peakon profiles along the diago- 
nals. However, the red (vertical) profiles show a symmetric progression, both 
upward and downward. This means the peakon strips are stretching as they 
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balloon out from the initial rightward motion. Being roughly circular at late 
times, this stretching of the strip length is comparable to the distance they 
propagate horizontally. Therefore, their stretching is an order 0(1) effect. 



Stability and an open problem. For a = a in Figure [19} the peakon 
curve segment is stable and it retains its integrity. The other cases, for a < a, 
are unstable and they break into narrower curved peakon segments each of 
width a, as the evolution tends toward the peakon profile e"'^'/". In fact, 
any smooth confined initial velocity distribution tends eventually to a set of 
peakon strips. These are contact curves, along which the discontinuity in 
slope moves with the fluid velocity. Hence, the momentum density tends to 
a set of moving curves on the plane, with each curve embedded in the flow. 
The latter claim is proved by the isospectral problem in Id, but for now 
is only an observation in 2d and 3d. The proof of this observed tendency 
remains an open problem from the analytical viewpoint. 



As a < cr decreases in Figures 20 22, the number of emerging contact 



curves increases. We emphasize that, at sufficiently late times, only these 
contact curves are observed emerging from a confined initial velocity distri- 
bution with width greater than a. However, the process occurs does require 
a certain amount of time to reach completion, as the peakons are successively 



formed at definite intervals. In the cases 20 and 21ofa = cr/2 and a = a/A 
at the times shown, some vestiges of the initial conditions still remain as 
ramps. As time progresses further, these ramps will eventually decay into 
a sequence of successively slower moving (lower amplitude) peakon contact 
curves. No such vestiges remain in the case when a = a. This behavior is 
in accord with the Id steepening lemma of [6], which states that an initial 
velocity profile possessing an inflection points of negative slope will develop 
a vertical slope in finite time. The formation of the vertical slope is part of 
the nonlinear steepening mechanism which creates the train of peakons from 
the ramp velocity configuration. 

Time reversal. Starting from the final velocity distribution for each value 
of a, we integrate back to the starting time in the EPDiff evolution numer- 
ics. This procedure tests the reversibility of the numerical algorithm. (The 
EPDiff equation itself is reversible.) Each case reverses accurately to its ini- 



tial condition, as is evident visually in Figure 23, and as measured in the 
L^, L^, and norms shown in Table [!} For the simple plate evolutions, 
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Figure 19: Evolution of the 2d "plate" initial velocity profile with a = a. 



without any of the head-on collisions of contact curves that will produce the 
more complicated flows seen below, this reversibility affirms the numerical 
scheme. It is not a severe test, however, in the sense that reversibility from 
the endpoint back to the beginning does not guarantee accuracy over the 
entire forward evolution. 

9.2 Parallel 



In Figures 24 27, two straight segments are initialized moving rightward. 
The one behind has twice the speed of the one ahead, and the two segments 
are offset in the vertical direction. The segments each break into curved 
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Figure 20: Evolution of the 2d "plate" initial velocity profile with a = a/2. 
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Figure 23: Reconstituted initial velocity profile after time reversal for the 
2d "plate" case with a = a (upper left), a = a/2 (upper right), a = a/A 
(lower left), a = (t/8 (lower right). This initial condition reconstitutes well 
for all values of a 
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strips of width a and these undergo overtaking coUisions. For a = a, the 
segments retain their integrity as they expand until the overtaking coUision 
occurs. Upon overtaking the slower segment, the faster segment transfers 
its momentum to the slower one ahead and a remarkable "reconnection" 
or "melding" of the segments occurs. This reconnection also shows rapid 
transverse stretching in which "hot spots" of velocity arise and then spread 
out along the wave front. Also, remarkably, the figures show a low amplitude 
"peakon wisp" connecting the endpoint of an earlier reconnection to a point 
(usually to a hotspot) on the leading peakon segment. This wisp apparently 
provides the "memory" of the previous reconnection, which is required for 
the evolution to remain reversible. Its reversibility affirms the nondissipative 
nature of the reconnection process. Moreover, this reconnection must be 
reversible, because the evolution is Hamiltonian in the continuum limit. 



The lower panels of Figures [24||27] for a = a show that the profiles remain 
exponential both before and after the overtaking collision. For decreasing 
values of a < cr, the evolution develops increasing complexity, with numerous 
overtaking collisions and corresponding reconnections. In each case, once sees 
the trailing memory wisps arising from these reconnections. These trailing 
wisps often connect to kinks at hotspots along the curve, indicating that 
their interaction is nontrivial, even though they have small amplitude. (This 
may also indicate that the hotspot is the source of the trailing wisp.) 



In the time-reversed runs shown in Figure 28, all of the overtaking col- 
lisions reconstituted their initial conditions accurately in the various norms 
shown in Table [Tj This indicates the accurate reversibility of the numerics 
for overtaking collisions. 

9.3 Skew 



Skew fiows in Figures 29 32 begin with two peakon segments of the same 
width, but oriented so that the one behind, which has twice the amplitude 
(speed) of the one in front, overtakes the one ahead by moving along the neg- 
ative diagonal. Again, one sees integrity of the a = a case and reconnection 
with hot spots and memory wisps trailing behind kinks in the main peakon 
segments after the overtaking peakon segment transfers its momentum to 
the one ahead. This locally Id soliton elastic-collision rule seems to explain 
the momentum transfer. Once again the lower panels show that the solution 
tends to peakon profiles in each direction. 
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Figure 28: Reconstituted initial velocity profile after time reversal for the 
2d "parallel" case with a = a (upper left), a = a/2 (upper right), a = a /A 
(lower left), a = cr/8 (lower right). 
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Memory wisps: another open problem. Figures 29 32 raises issues 
(such as the production of memory wisps) for the 2d peakon segment in- 
teractions which could lead to new research well beyond the scope of the 
present work. As before, in the cases of a < a, for skew flows it seems that 
each memory wisp is attached to a hot spot along a major peakon segment. 
The memory wisps, by the way, are very low in amplitude, but they seem to 
also possess the peakon exponential transverse profile. Hence, the stretching 
motion along the wisp must considerably faster than across it. This may be 
because the hot spots keep contributing to the wisps trailing behind them. 
Thus, the wisps may be seen as trailing residuals from the hot spots. The 
memory wisp feature of the reconnection remains to be explained in more 
detail, both numerically and analytically. 



For small a < a, the central regions of the skew flows in Figures 30 



32 become very intricate (partially mixed). Hence, one could expect that 
its reversibility is compromised. The case a = a/8 is mostly reconstituted 
after the time reversal, but some small-amplitude, high-frequency errors do 



remain, as observed in Figure 33 For other values of alpha {a = a , a/2 , a /A) 



the initial conditions reconstitute quite well. 



9.4 Wedge 



The wedge flows in Figures |34H37| are variants of skew flows in Figures [29p2 
showing two plates colliding along opposite diagonals in the plane, with re- 
flection symmetry about the horizontal axis. The wedge flows are convergent, 
and therefore they have some head-on features that emerge on the left-hand 
side of the collisions. They also show considerable acceleration along the 
midline, in forming jets moving along the horizontal axis in both directions. 
This jet formation is due to convergence of momentum which continues to 
build up after the initial collision. The multiple wedge collisions occurring 
for values of a < a show successive strong accelerations due to convergence. 
They also show enhanced stretching of the main peakon segments. These 
collisions also produce complex patterns of small-amplitude peakon wisps, 
trailing from hotspots at kinks along the main peakon curves. 



Again, the lower panels in Figures 29 32 show primarily peakon profiles 



and the emergence of peakons from ramp-and-cliff formations with inflection 
points of negative slope, in agreement with the steepening lemma for CH in 
Id. 

The head-on features of the wedge collisions cause noticeable errors in 
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Figure 29: Evolution of the 2d "skew" initial velocity profile with a = a. 
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Figure 30: Evolution of the 2d "skew" initial velocity profile with a = a/2. 
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Figure 31: Evolution of the 2d "skew" initial velocity profile with a = a/4. 
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Figure 32: Evolution of the 2d "skew" initial velocity profile with a = a/8. 
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Figure 33: Reconstituted initial velocity profile after time reversal for the 
2d "skew" case with a = a (upper left), a = a/2 (upper right), a = a/A 
(lower left), a = a/S (lower right). 
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Figure 34: Evolution of the 2d "wedge" initial velocity profile with a = a. 



reversibility, of order 10 to 20 percent for a = cr/8, as observed in Figure 



38| However, these errors decrease for larger alpha. As we will see in other 
plots, errors in reversibility tend to arise for smaller values of a < cr whenever 
head-on collisions occur. 



9.5 Head-on 



The head-on collisions of two offset peakon segments in Figures 39 42 show 
great complexity. Some of this complexity is due to the process of annihila- 
tion and recreation known to occur in the purely Id antisymmetric head-on 
collisions of a peakon with its reflection, the antipeakon. At the moment 
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Figure 36: Evolution of the 2d "wedge" initial velocity profile with a = a/4. 
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Figure 37: Evolution of the 2d "wedge" initial velocity profile with a = a/8. 
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Figure 38: Reconstituted initial velocity profile after time reversal for the 
2d "wedge" case with a = a (upper left), a = a/2 (upper right), a = a/A 
(lower left), a = cr/8 (lower right). 
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a head-on collision occurs, one sees the blue profiles in the lower panels of 



Figures [39H42] becoming very steep, exactly as seen in Id peakon-antipeakon 
collisions. Thus, the locally Id collision rules carry over to the 2d case for 
head-on collisions. 



The case a = a in Figure 39 shows the annihilation and recreation ex- 
pected from the Id rules. Thereafter, a new feature emerges as the recreated 
segments disconnect, then eventually reconnect again with segment elements 
that did not participate in the head-on part of the collision. The first discon- 
nection occurs in the presence of rapid rotation or circulation of the peakon 
segments, which causes extreme stretching before the disconnection. The 
later reconnection of these segments occurs less violently. 

Once again, peakon profiles are ubiquitous in the lower panels of Figures 

Time reversal of the head-on collisions shows, for a = a, an essentially 



complete restoration of the initial condition; see Figure 43 However, the 
cases for smaller alpha [a = a/4 ,a/8) show breakup and head-on collisions 
occurring simultaneously. These simultaneous processes produce complex 
mixed states which tend to reverse less accurately, as one might expect, 
because of the plethora of peaked excitations at small scales. 

9.6 Star 



Star flows in Figures [44]|47] form a variant of wedge flows with fivefold sym- 
metry, instead of simple reflection antisymmetry. The mutual rotation of the 
overtaking-collision evolution in each case preserves the fivefold symmetry 



well, and is seen in Figure 48 to be largely reversible. Figures 44 47 each 
show many reconnections (mergers), until eventually one peakon filament 
ring surrounds all the others. Again, peakon segments are the ubiquitous 
feature of the solution. If the evolution were allowed to proceed further, 
reconnections would tend to produce additional concentric rings of peakon 
filaments. 



9.7 Rotate 



Figures 49 52 show how angular (azimuthal) motion couples to radial motion 
in the plane. Each evolution starts with a circularly symmetric velocity 
distribution in a Gaussian (not peakon) ring, or annulus, of width alpha, 
which is initially rotating rigidly at constant angular velocity. The angular 
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Figure 43: Reconstituted initial velocity profile after time reversal for the 
2d "head-on" case with a = a (upper left), a = a/2 (upper right), a = a /A 
(lower left), a = cr/8 (lower right). 
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Figure 45: Evolution of the 2d "star" initial velocity profile with a = cr/2. 
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Figure 46: Evolution of the 2d "star" initial velocity profile with a = cr/4. 
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Figure 48: Reconstituted initial velocity profile after time reversal for the 
2d "star" case with a = a (upper left), a = a/2 (upper right), a = a/A 
(lower left), a = cr/8 (lower right). 
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motion couples to the radial motion by producing a sequence of outward 
and inward propagating circular peakons. These circular peakons also rotate 
clockwise (in the same sense as the initial condition) because of conservation 
of angular momentum. See Holm et al. [25|[26] for more details about circular 
peakons. 



For a < a, the first inward propagating circular peakon in Figures 50 52 



collapses to the center, then reflects outward and overwhelms the formation 
of any subsequent inward moving circular peakons. For a = cr, no inward 
moving peakons form during the time of the simulation. 

Upon reflection from the center, the circular peakon is influenced by the 
finite Cartesian mesh, especially for smaller alpha (narrower peakons). The 
refiection interaction with the mesh distorts the outward propagating solu- 
tion, as seen in the animations, with greater distortion for smaller alpha. This 
is a severe test of the numerics near the center of symmetry. Nonetheless, 
time reversal from the final velocity profile reconstitutes the initial condition 



to within less than one percent, as seen in Figure 53 This illustrates that 



accurate time reversal to the initial condition does not guarantee accuracy of 
the solution during the forward evolution. However, it does show that even 
in a somewhat distorted solution, dissipation plays little role in the numerical 
simulation. 



9.8 Right 



In the simulations shown in Figures 54 57, the velocity distribution is ini 



tially a Gaussian ring in magnitude, uniformly pointed rightward along the 
X axis. The right outer side of the ring produces diverging peakon contact 
curves, which slow as they propagate outward. The left inner side of the 
ring, however, produces converging peakon contact curves, which accelerate 
as they converge, undergo a strong interaction along the axis, then break 
again into contact curves still moving rightward, approaching the previous 
divergent peakon curves and colliding with them from behind. These over- 
taking collisions impart momentum but they do not produce reconnections. 

After the collisions, a complex fiow remains near the axis, in which one 
also sees hot spots at kinks in the contact curves, with trailing memory wisps 
behind them. 

The lower panels show peakon profiles with high wavenumber oscillations 
(possibly noise) in the complex fiow region remaining behind near the axis. 
Except for the smallest case of a = cr/8, all of the time-reversed runs re- 
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Figure 49: Evolution of the 2d "rotate" initial velocity profile with a = a. 
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Figure 50: Evolution of the 2d "rotate" initial velocity profile with a = a/2. 
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Figure 52: Evolution of the 2d "rotate" initial velocity profile with a = a/8. 



D. D. Holm & M. F. Staley 



Singular Wave Fronts 



88 




Figure 53: Reconstituted initial velocity profile after time reversal for the 
2d "rotate" case with a = a (upper left), a = a/2 (upper right), a = a/A 
(lower left), a = a/S (lower right). 
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Figure 54: Evolution of the 2d "right" initial velocity profile with a = a. 



assemble into the Gaussian ring without significant distortion, as diagnosed 



by the black profile in the lower panel of the time reversed runs in Figure 58 



9.9 Inout 



In Figures 59 62 we start with an initial Gaussian ring of width alpha in 
speed, with an angular distribution of — (sin6', cos 6^) exp(— (r — ro)^/cr^) for 
the direction of the velocity. Consequently, the motion is inward along the 
positive diagonal and outward along the negative diagonal. The outward 
motion breaks into a sequence of curved peakon segments of width alpha, as 
usual. The inward motion also produces peakon segments, which however 
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Figure 55: Evolution of the 2d "right" initial velocity profile with a = a/2. 
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Figure 58: Reconstituted initial velocity profile after time reversal for the 
2d "right" case with a = a (upper left), a = a/2 (upper right), a = a/A 
(lower left), a = a/S (lower right). 
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Figure 59: Evolution of the 2d "inout" initial velocity profile with a = a. 



undergo head-on collisions, so that they annihilate, recreate and re-emerge 
moving along the positive diagonal. 



The blue profiles in the lower panels of Figures [59H62] show the breakup 
of the outward motion into peakon profiles. The green profiles on the lower 
panel show the head-on collisions of the peakon profiles of the inward motion. 
The inward moving head-on collisions leave a residue of complex flow. 

Time reversal in this case shows severe distortion of the initial condition 
for the smaller values a = a/8 and a = (t/4, but not for the larger values 
a = a/2 and a = a. See Figure 63 
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Figure 60: Evolution of the 2d "inout" initial velocity profile with a = (t/2. 
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Figure 61: Evolution of the 2d "inout" initial velocity profile with a = a/4. 
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Figure 63: Reconstituted initial velocity profile after time reversal for the 
2d "inout" case with a = a (upper left), a = a/2 (upper right), a = a/i 
(lower left), a = a/8 (lower right). 
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9.10 Time reversal 

Table [T] summarizes how well the initial velocity profiles reconstitute upon 
reversing the EPDiff evolution from the final states. Numbers in each table 
entry are the norm, norm, and max norm of the difference between the 
initial velocity profile and its reconstituted (time reversed) value. The time 
reversed values are typically accurate to within one percent or less except in 
cases where head-on collisions occur. See text for further discussion. 

10 Numerical results for EPDifF in 3d 

Singular solutions of EPDiff in 3d. In this section, we extend our nu- 
merical solutions of EPDiff to 3d and discover that the codimension-one 
singular solution behavior of EPDiff persists. Namely, the singular solution 
behavior of EPDiff extends from curves in 2d to surfaces in 3d. To our 
knowledge, this is the first example of numerical simulations of the nonlinear 
interactions of contact discontinuities for fluid velocity in 3d. Of course, the 
singular surface solutions of EPDiff in 3d have no interpretation in the the- 
ory of internal waves. However, codimension-one singular solutions of EPDiff 
in 3d may have applications elsewhere, for example in the physics of liquid 
crystals in the inertial regime, as discussed for Id evolution by Hunter and 
Saxton [30]. The 3d evolution of contact surfaces may also be regarded as 
"growth" from the viewpoint of imaging science, or evolution of shapes in 
3d. The present work focuses on the role of convergence and momentum 
exchange in the evolution of contact surfaces. 

Upcoming sections describe 3d simulations of evolution under EPDiff for 
each of twelve initial velocity distributions. For each simulation we have 
figures for a = a and a = a/2. Due to memory and computational limi- 
tations, our 3d simulations were run at one-fourth the resolution of our 2d 
simulations (256'^ instead of 1024^), and this lower resolution precluded us 
from performing accurate computations at the smaller values a = a/A and 
a = a/8 that were simulated in 2d. 

Each figure contains six frames, showing the initial magnitude (speed) 
of velocity, |u|, in the upper left frame, followed by plots of |u| at future 
times, reading across and then down. For each simulation, the domain is 
[— 1, 1] X [—1, 1] X [—1, 1], and x is toward the right, y toward the back, and z 
toward the top. Individual plots contain colored, partially transparent level 
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Table 1: L^, L'^, and max norm of the difference between the initial velocity 
profile and its reconstituted (time reversed) value for each 2d simulation. 



Simulation 


a = a 


a = a/2 


a = a/4: 


a = a/8 


plate 


3 . 62e-05 
0.000142 
0.00296 


7 . 03e-05 
0.000199 
. 00279 


. 000213 
0.000624 
0.00371 


. 000577 

0.00207 

0.0147 


parallel 


7 . ooe-Oo 
0.000204 
0.00301 


8 . 69e-05 
0.000217 
0.00283 


0.00024 

0.000616 
0.00358 


0.000844 

0.00236 

0.0161 


skew 


. 00134 
0.0029 

0.0199 


. 000498 
0.00111 

0.0102 


. 000719 
0.00139 

0.00727 


. 00711 
0.0177 

0. 148 


wedge 


. 000532 
0.000705 
0.00395 


. 00061 
0.00108 
0.00539 


. 00451 
0.0123 
0. 148 


. 0567 
0. Ill 
0.671 


head-on 


. 000179 
. 000332 
. 00335 


. 00636 

0.0165 

0.172 


0.0274 
0.0622 
0.43 


0.0538 

0.111 

0.665 


star 


0.000207 
0.000419 
0.00531 


0.00023 

0.000442 

0.00476 


0.000604 
0.000851 
0.00515 


0.00152 
0.00237 
0.00883 


rotate 


7.21e-05 
8.17e-05 

0.000301 


2 . 54e-05 
3.95e-05 

0.000139 


4.81e-05 
8.23e-05 

0.000339 


0.000293 

0.00059 

0.0027 


right 


0.00021 
0.000372 
. 00545 


0.000292 
0.000615 
. 00422 


0.000583 

0.00137 

0.0106 


0.0029 

0.00696 

0.0892 


inout 


0.000219 

0.00048 

0.00621 


. 000245 

0.000431 
0.00183 


0.0024 

0.0103 
0.137 


0.0348 

0.087 

0.594 
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Figure 64: Locations of 2d slices of |u| shown in the left, back, and bottom 
panels in the upcoming 3d figures. 



surfaces of |u| at 35% (purple), 25% (red), 15% (green), and 5% (blue) of its 
maximum value over all frames in the figure. The number of level surfaces, 
and the colors, amount of transparency, and fractions of |u| represented by 
each level surface, where chosen after considerable experimentation so that 
they reveal maximal information with minimal clutter. 

The left, back, and bottom panels of each frame show color contour plots 
of 2d slices through the three dimensional data at x = 0, y = 0, and z = 0, 



respectively, as sketched in Figure 64 The same color scheme that was used 



for the 2d simulations, as described in Section [9] and illustrated in Figure 



17 is used here for the 2d slices of the 3d simulations. Note that the colors 
used for the 2d slices do not correspond to the colors used for the 2d level 
surfaces. For example, red in the 2d slices does not represent the same value 
of lul as the red level surface. 



10.1 Plate 



The analog of a segment in 2d is a disc, or plate, in 3d. Figures 65 and 66 
show the evolutions at various values of alpha in 3d from a confined initial 
velocity distribution in the shape of a plate in speed and moving rightward 
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(in the positive x direction). The initial plate distribution is chosen so it 
falls off exponentially in the normal x direction and it falls off at the edges 
(or rim) of the plate as a Gaussian, as we did in 2d at the endpoints of the 
segments. Because the distribution moves with flow, it expands rightward 
and outward. 

The bottom panels in Figures 65 and 66 show the propagation in a hor- 
izontal section (at the midplane) x = of the cube of size 2. This is an 
invariant plane, by symmetry, for this initial value problem. The invariant 
midplane section essentially reproduces the 2d evolution of the peakon seg- 
ments; note the similarity in patterns at late time on the bottom panel in 3d 
and the corresponding interactions in 2d. 

The bottom panel in Figure 66 with a = a/2 shows the plate expanding 
and decomposing into several peakon contact surfaces. The left panel shows 
a vertical section at x = 0. The initial circular disc propagates rightward 
through this vertical section and maintains its circular symmetry as it ex- 
pands. The back panel shows a vertical section aX y = 0. By symmetry, the 
results on the vertical section shown on the back panel in this case are the 
same as those on the bottom panel {z = horizontal section). 

The disc shaped velocity distribution balloons out rightward and expands 
cylindrically into one peakon contact surface in Figure 65 for a = a, and into 
several peakon contact surfaces of width a = cr/2 in Figure 66 Again, this is 
consistent with the similar evolution restricted to lower spatial dimensions. 



10.2 Parallel 



Figures 67 and 68 each show two plate-like velocity distributions of the same 
diameter, initialized moving rightward so that the one behind, with twice 
the speed, will overtake the one ahead. The discs are initially offset at 45 
degrees in the y-z plane; so this evolution has no reflection symmetry about 
any of the midplanes. 

Both plates balloon outward and decompose into peakon contact sur- 
faces, and the overtaking collision shows two reconnections of the surfaces. 
One reconnection occurs at first contact and the other evolves late in the 
simulation, as seen on both the bottom and back panels. The left panel 
shows a transverse vertical section of the first reconnection. As the two discs 
propagate past the plane x = 0, their peakon contact surfaces develop into a 
perimeter that surrounds the interior interactions. 
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Figure 65: Evolution of the 3d "plate" initial velocity profile with a = a. 




Figure 66: Evolution of the 3d "plate" initial velocity profile with a = a /2. 
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10.3 Skew 



The skew overtaking collision in Figures |69] and [70] show both the first and 
second reconnections as the two initial plate distributions collide and interact. 
Again, the sections shown on the bottom and back panels are very reminiscent 
of the corresponding skew collisions of peakon segments in 2d. The memory 
wisps that appear so clearly in the 2d skew overtaking collisions are less 
distinct here, likely because of the coarser resolution (256^, versus 1024^ in 
2d). 

10.4 Wedge 



In Figures [71] and [72] two plate distributions of velocity approach the z = 
horizontal midplane from above and below, each at 45 degrees. Their 
"wedge" collision occurs at the horizontal midplane and has reflection sym- 
metry about the xz vertical midplane, at y = 0. 

The evolution shown on the back panels in these figures is similar to the 
2d wedge collision of peakon segments. The bottom panel shows the creation 
of peakon contact surfaces which are reminiscent of the peakon contact seg- 
ments in the 2d "right" figure. Namely, peakon contact surfaces emerge from 
the collision and reconnect to encircle the collision region. Each successive 
peakon contact surface (whose intersection with the horizontal midplane is 
shown as a curve in the bottom panel) reconnects at the rear and produces a 
wavefront that encircles the collision region. The wavefront motion is less in- 
tense in the rear because it is moving slower at the rear than in the front. We 
see one major collision at the midplane followed by an emission of rightward 
moving peakon contact surfaces and their later reconnections in the rear. At 
later times, the hot spots show pronounced trailing memory wisps. 

10.5 Head-on 



Figures 73 and 74 show the head-on collision of two identical plate-like dis- 
tributions offset by the same distance at 45 degrees when projected into the 
yz plane. Hence, the evolution in the sections shown in the back and bot- 
tom panels are identical, and all three panels have reflection and rotation 
symmetries. 



The back and bottom panels of Figures 73 and 74 have symmetry under 



reflection about one diagonal and rotation by vr about the other. Likewise, 




Figure 69: Evolution of the 3d "skew" initial velocity profile with a 



= a. 
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Figure 70: Evolution of the 3d "skew" initial velocity profile with a = a/2. 




Figure 71: Evolution of the 3d "wedge" initial velocity profile with a = a. 




= (t/2. 
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the left panels are reflection symmetric about both diagonals. These sym- 
metries help diagnose the complex evolution occurring as the two peakon 
contact segments balloon outward toward each other and collide head-on. 
As expected from earlier animations, the segments annihilate and re-emerge, 
leaving behind a complex residual flow where the head-on collision occurred. 

The left panels of Figures [73] and 74 show that the reconnection after the 
collision occurs results in a peakon contact segment surrounding the central 
region. 



10.6 Rotate 



For the simulations shown in Figures 75 and 76 , an initially spherical Gaus- 



sian shell is rigidly rotating about its vertical axis, like a planet. This angular 
motion couples to radial motion with cylindrical symmetry about the axis 
of rotation. The evolution remains cylindrically symmetric about this axis 



and shows convergence to it, as seen in the bottom panels of Figures 75 and 



76 The left and back panels show identical emergence of peakon contact 
segments. This evolution shows the coupling between the angular and radial 
motion in the cylindrically symmetric case. The outward velocities expand 
essentially like an oblate sphere, and may become less oblate with time. The 
midplane z = shown on the bottom panel is an invariant section, because of 
up-down symmetry, and it shows the cylindrical convergence and expansion 
driven by the initial angular rotation. 



10.7 Right 

In Figures [77] and [78} a Gaussian shell in speed is initially moving rightward. 
The peakon contact surfaces emerging on its outer right side are diverging, 
while those emerging on its inner left side are converging. The acceleration 
due to this convergence leads to an overtaking collision that imparts right- 
ward momentum to the diverging peakon contact segments. By cylindrical 
symmetry about the x axis, the bottom and back panels show the same mo- 
tion. The left panels show how this cylindrically symmetric motion moves 
through the vertical midplane at x = 0. The midplane z = is invariant and 
the motion in this plane mimics the motion observed in the corresponding 
2d problem. 

This shared behavior is one of the main conclusions of this paper: 3d 
numerical results have planar slices which show the corresponding 2d mo- 
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Figure 73: Evolution of the 3d "head-on" initial velocity profile with 
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Figure 75: Evolution of the 3d "rotate" initial velocity profile with a = a. 
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Figure 76: Evolution of the 3d "rotate" initial velocity profile with a = a/2. 
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mentum transfer behavior, and 2d numerical results have linear slices which 
show the corresponding Id momentum transfer behavior. This reduction 
principle allows the complex interactions of contact wave surfaces to be ana- 
lyzed as elastic collisions showing local momentum transfer. This is a feature 
of the collective behavior of singular solutions, because such solutions possess 
no internal degrees of freedom. 



10.8 Inout 



In Figures [79] and |80} a 3d Gaussian shell of speed is initially given the 
2d "inout" velocity distribution, — (sin6', cos6') exp(— (r — roY/cr'^), on each 
horizontal plane, weighted by the cosine of the polar angle, as in the rigid 
rotation simulation. Consequently, the bottom panel shows the same "in- 
out" 2d motion as before, now at one-fourth the resolution and only for 
a = a/2 and a = a. The left and back panels show identical motion, 
consisting of head-on collisions in the center and a wedge-like collision at the 
top and bottom. The bottom panel for a = a/2 seems to agree better with 
the corresponding 2d case, than for a = a, which has less of a collision that 
in the 2d case. 



10.9 Wheel 



In Figures [81] and [82] the speed is initially distributed as a Gaussian within 
a toroidal annulus rotating rigidly about its axis of rotational symmetry, ori- 
ented along the x axis. This "wheel" initial distribution has reflection sym- 
metry about its axis of rotational symmetry, oriented along the midplane 
a; = 0. The motion on this invariant symmetry plane x = shows circularly 
symmetric collapse and radial expansion, seen in the left panel. (Some dis- 
tortion is seen upon reflection from the x axis of cylindrical symmetry, again 
due to mesh effects.) 

The identical back and bottom panels show a head-on collision occurring 
on the horizontal midplane, followed by re-emergence in 3d for a = a/2 and 
the development of a pair of strong hot spots (actually a funnel shape in 3d) 
followed by emergence of a perimeter of peakon contact segments surrounding 
the interior region. 
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Figure 77: Evolution of the 3d "right" initial velocity profile with a = a. 




Figure 78: Evolution of the 3d "right" initial velocity profile with a = a/2. 
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Figure 79: Evolution of the 3d "inout" initial velocity profile with a = a. 
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Figure 80: Evolution of the 3d "inout" initial velocity profile with a = a/2. 




Figure 81: Evolution of the 3d "wheel" initial velocity profile with a = a. 
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Figure 82: Evolution of the 3d "wheel" initial velocity profile with a = a/2. 
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10.10 Wheels 



Figures [83] and [84] show simulations of the interactions of two coaxial tori 
(wheels) of velocity, both rotating about the x axis, but offset from each 
other in the x direction. The cylindrically symmetric evolution from this 
initial state shows axial convergence, radial expansion and reconnection of 
peakon contact surfaces, via a series of head-on collisions. The exchange 
of momentum in these interactions is especially dramatic. Eventually, an 
outward-expanding peakon contact surface emerges and surrounds the entire 
interaction region. 

10.11 Torus 



In Figures 85 and 86 , a torus-shaped Gaussian velocity distribution (wheel) 
is initially moving uniformly rightward in the direction of its symmetry axis. 
Cylindrically symmetric rightward moving fronts form, then expand, collide, 
and reconnect at the axis of symmetry. This is a "cylindrically symmetric 
wedge" collision that funnels into the axis, then forms jets that accelerate in 
both forward and backward directions along the 

The three dimensional images of peakon contact surfaces are particularly 



vivid in the animations to which Figures 85 and 86 belong. Because of the 
cylindrical symmetry, the back and bottom panels show the same motion 
in different perspectives. The left panel also shows the collapse to the axis 
of symmetry. Because of the curvature of the peakon contact surfaces, this 
cylindrical collision imparts axial momentum both forward and backward, 
which is especially clear in the case a = a/2. 

10.12 Tori 



In Figures 87 and 88 two linked tori of speed are started along their axes of 
rotational symmetry in orthogonal directions, one rightward and one upward. 
They undergo a series of wedge-like collisions leading to many reconnections. 
They also undergo overlapping collisions that impart momentum, but do not 
reconnect the peakon contact segments. Eventually, a single outward moving 
peakon contact segment will surround the interaction region. 



D. D. Holm & M. F. Staley 



Singular Wave Fronts 




Figure 83: Evolution of the 3d "wheels" initial velocity profile with 




Figure 84: Evolution of the 3d "wheels" initial velocity profile with a = a /2. 




Figure 85: Evolution of the 3d "torus" initial velocity profile with a = a. 




Figure 86: Evolution of the 3d "torus" initial velocity profile with a = a/2. 
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Figure 87: Evolution of the 3d "tori" initial velocity profile with a = a. 
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Figure 88: Evolution of the 3d "tori" initial velocity profile with a = a/2. 
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10.13 Time reversal 

Table [T] summarizes how well the initial velocity profiles reconstitute upon 
reversing the EPDiff evolution from the final states. Numbers in each table 
entry are the norm, norm, and max norm of the difference between 
the initial velocity profile and its reconstituted (time reversed) value. We 
do not show contour plots of the time reversed runs in 3d because they 
are all visually indistinguishable at the larger values of a chosen for our 3d 
simulations. The time reversed values are typically accurate to within one 
percent or less. 

11 Conclusions, future directions, and out- 
standing problems 

By a sequence of approximations applied to the Euler-Poincare variational 
principle for the multilayer columnar motion (MLCM) of an incompressible 
fluid, we derived a hierarchy of EP equations. Several new MLCM equa- 
tions belong to this hierarchy, as well as the standard Boussinesq equations 
for shallow water waves, which are recovered upon specialization to weak 
nonlinearity. 

The EPDiff equation Q was derived in the limiting case of a single layer 
undergoing strongly nonlinear motion in the absence of linear dispersion. In 
Id, the EPDiff equation restricts to the dispersionless case of the CH equation 
for nonlinear shallow water waves. The dispersionless CH equation possesses 
singular soliton solutions whose velocity possesses a sharp peak (jump in 
slope) moving at a speed equal to its height. The corresponding momentum 
for a train of such peakons is a set of delta functions at the locations of 
the peaks in velocity. Thus, this momentum density is distributed as points 
on the line which evolve under the action (by Ad*) of the diffeomorphisms 
(smooth invertible maps). 

A geometrical version of the soliton paradigm. The local description 
of the Ad* action of the smooth invertible maps is the EPDiff ad* equation 
([T]), which holds in any number of dimensions. Hence, EPDiff allows compar- 
ison of the behavior of its singular solutions in Id, 2d and 3d. Numerically, 
the singular solutions of EPDiff of codimension one (points on the line, curves 
on the plane, surfaces in a volume) are found to be stable. Moreover, they 
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Table 2: L^, L^, and max norm of the difference between the initial velocity 
profile and its reconstituted (time reversed) value for each 3d simulation. 



Simulation 


a — a 


a = a/2 


plate 


0.000105 
0.000595 
0.0109 


9.26e-05 
0.000553 
0.01 


parallel 


0.000165 
. 000686 
0.0109 


0.000146 
. 00064 
0.01 


skew 


9.81e-05 
. 000435 
. 00793 


8.68e-05 
0.000411 
. 00753 


wedge 


0.00029 

0.000513 
0.00545 


0.000135 

0.000398 
0.00522 


head-on 


0.000193 
0.000772 
. 00885 


0.000178 
0.000738 
. 00829 


rotate 


9.36e-05 

0.000151 
0.000625 


6.13e-05 

0.000102 
0.000433 


right 


6.08e-05 
0.000128 
0.00131 


4.01e-05 
8.88e-05 
0.00103 


inout 


9.4e-05 

0.000151 

0.00063 


6.19e-05 
0.000103 
0.000436 


wheel 


5 . 87e-05 
0.000176 

0.0019 


3.98e-05 
0.000124 
0.00138 


wheels 


0.000123 
0.000259 
0.00193 


8.51e-05 
0.000184 
0.0014 


torus 


7.92e-05 
0.000109 
0.000523 


5.45e-05 
7.87e-05 
0.000391 


tori 


0.00011 

0.000248 

0.0019 


7.48e-05 
0.000173 
0.00137 
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are found to dominate the intial value problem, essentially by following the 
soliton paradigm. Namely, the "singular solution content" of a given initially 
continuous distribution of velocity emerges under the evolution of EPDiff and 
retains its integrity under collision interactions. In Id, these singular solu- 
tions are true solitons (the peakons) for the dispersionless CH shallow water 
equation. In this case, the inverse scattering transform for CH determines 
its soliton content for an arbitrary initial condition and also gives the soliton 
collision rules. In 2d and 3d, the theory for determining how these singular 
solutions will emerge and how they will survive their collisions has not yet 
been developed. This is the major open problem for EPDiff. One clue for 
approaching it is to keep in mind that the singular solutions in ^ form an 
invariant manifold for the EPDiff equations. 

Our numerical findings. We studied the initial value problem for EPDiff 
in various scenarios in 2d and 3d. We found that the key solution behavior 
of the IVP for EPDiff is the breakup of an initially smooth confined veloc- 
ity distribution into singular solutions supported on codimension-one delta 
function densities moving with the velocity of the flow. We extended the so- 
lution ansatz for peakon momentum density supported on points on the line 
to the case of singular momentum density supported on smoothly embedded 
sets, e.g. curves in the plane, or surfaces in three dimensional space. We 
showed that this singular solution ansatz for EPDiff reduces to Hamilton's 
canonical equations for the vector parameters defining these surfaces. The 
underlying geometrical reason why this reduction occurs was explained in 
Holm & Marsden [22] by recognizing that the singular solution ansatz ^ is 
a momentum map for the (right) action of diffeomorphisms on distributions 
defined as smoothly embedded subspaces of a manifold. Thus, the singular 
solutions evolve by Ad* action on embedded subspaces along a curve g{t) in 
the diffeomorphisms which is a geodesic path. As we explained, such a curve 
is a geodesic if and only if the corresponding momentum satisfies the EPDiff 
equation ([T|. 

Remarkably, our numerical results showed that only the codimension- 
one singular solutions emerge in the IVP. Being defined on delta functions, 
these solutions have no internal degrees of freedom. Consequently, their local 
interactions may be characterized as elastic collisions of contact surfaces in 
which momentum is exchanged. Across these contact surfaces, the slope of 
the velocity has a jump which moves with the flow. The collision rules for 
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these interactions in 2d and 3d may be built up from the sohton colhsion 
rules in Id. That is, a linear section transverse to a 2d solution shows Id 
elastic momentum exchange behavior, and a planar section transverse to a 
3d solution shows the corresponding 2d behavior. This reduction to lower 
dimensional behavior holds especially well on reflection-invariant sections. 
For example, the midplane z = is invariant in Figure 78 and the motion 
in this plane projected on the bottom panel of Figure 78 mimics the motion 



observed in the corresponding 2d problem in Figure 55 



Thus, 3d numerical results have planar slices which show the correspond- 
ing 2d momentum transfer behavior. Likewise, 2d numerical results have 
linear slices which show the corresponding Id momentum transfer behavior. 
Such a reduction principle allows the complex interactions of contact wave 
surfaces to be analyzed as elastic collisions showing local momentum transfer. 
This principle for collective behavior based on simple momentum exchange in 
collision interactions arises for the singular solutions, because such solutions 
possess no internal degrees of freedom. 

Two new numerical features emerge in 2d and 3d. These are the re- 
connections of the wavefronts, which occur due to momentum exchange in 
these nonlinear collisions, and the remarkable "memory wisps" that arise to 
guarantee reversibility of those collisions. (These memory wisps are a bit 
reminiscent of the emission of neutrinos, which preserve detailed balance in 
beta decay.) The memory wisp feature of the reconnections remains to be 
explained in more detail, both numerically and analytically. 



EPDiff applications: solitons, turbulence and medical images. 

From the viewpoint of nonlinear PDE analysis, the EPDiff equation, be- 
ing nonlinear and nonlocal, escapes classification. Its nonlocality requires 
solving an elliptic problem for determining its velocity from its momentum 
at each time step. It's worth repeating that the EPDiff nonlinearity is ge- 
ometric, because this is the key to understanding its motion. Namely, it is 
reversible geodesic motion in the diffeomorphisms acting on smoothly embed- 
ded subspaces. Physically, the singular solutions are contact surfaces (jumps 
in the velocity derivative that move with the flow). The corresponding EPDiff 
equation on the volume-preserving diffeomorphisms is the Lagrangian Aver- 
aged Euler (LAE-a) equation derived first in [23], which was first derived 
as a geometrical extension of the CH equation. Upon adding Navier-Stokes 
viscosity, this became the LANS-a model for incompressible turbulence in 
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[2] . For a review of the properties of LANS-a solutions and their relation to 
Navier-Stokes analysis, see ^J] and [S]. 

By a remarkable coincidence, L. Younes derived the EPDiff equation ([T]) 
as the evolution for the template matching approach in medical imaging in 
|52j . See also D. Mumford's discussion of the same problem in |1H]. Thus, the 
singular solutions we discuss here, and their momentum exchange paradigm, 
should be expected to develop increasing interest in medical imaging science. 
Medical imaging must solve an optimization problem rather than the IVP. 
However, our solution of the IVP may help guide intuition in medical imaging, 
especially the idea of encoding information in an image by its momentum, 
rather than just by the positions of its outlines and their associated intensi- 
ties. Thus, EPDiff represents a crossroads of endeavor in mathematics where 
methods of fluid dynamics and imaging science may transfer technologies. 
See [27] for more discussions of this new paradigm in image processing. 

Many open questions remain. Many open problems and other future 
applications remain for the EPDiff equation. For example, its analysis re- 
quires additional developments of PDE methods. In particular, while its 
smooth solutions satisfy a local existence theorem analogous to the Ebin- 
Marsden theorem for the Euler fluid equations, its IVP inevitably develops 
singular solutions. The implications of this observation are mentioned in 
[22] as perhaps indicating an incompleteness of the geodesic flows on the 
diffeomorphisms, which opens many future opportunities for analysis of the 
emergence of measure values solutions from smooth initial conditions in non- 
linear nonlocal PDEs. 

In addition, many open questions remain for the practical problem of in- 
ternal wavefronts emerging in shallow water dynamics. For example, there re- 
mains the issues of boundary and topography interactions, including diffrac- 
tion and refraction. Moreover, for the applications in turbulence modeling 
there remains a variety of open questions about singular vortex interactions, 
which are responsible for the famous vortex stretching that drives the cas- 
cade of energy and vorticity in turbulence. A host of other problems also 
remains for the applications of EPDiff in medical imaging, particularly for 
the statistical treatment of the information encoded in the linear space of 
their momenta. Finally, the development of numerical approaches that are 
fully capable of tracking the singular solutions of EPDiff, perhaps by using 
geometrical methods which incorporate discrete exterior calculus and varia- 
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tional multisymplectic integration methods. See [12] and [SB] for descriptions 
of these promising methods, which seem to he on the horizon for the next 
apphcations of singular solutions. 

Acknowledgements 

We are grateful to many friends and colleagues for their enormous help and 
encouragement during the course of this work. These include R. Lowrie 
and B. A. Wingate, whose encouragement inspired the beginning of this en- 
deavor in the context of Lagrangian averaged models of turbulence. We are 
also grateful to R. Camassa, A. Hirani, M. Leok, S. T. Li, J. E. Marsden, 
M. I. Miller, T. Ratnanather, and L. Younes for helpful discussions and chal- 
lenging observations. Finally, we thank the Los Alamos Turbulence Working 
Group (TWG) and the members of the TWO external advisory committee 
for continuing discussions, challenges, and encouragement. 

This work was supported by the United States Department of Energy 
under contracts W-7405-ENG-36 and the Applied Mathematical Sciences 
Program KC-07-01-0L We are also grateful for partial funding from Applied 
Mathematical Sciences Division, Office of Advanced Scientific Computing 
Research, DOE Office of Science. 

References 

[1] This process may be extended to higher moments, at the cost of intro- 
ducing additional variables, which we decline. 

[2] Note: The top and bottom boundary conditions influence stability char- 
acteristics. 

[3] Contributions arising when integrating by parts vanish, upon requiring 
the usual condition that the fluid velocities and their variations are tan- 
gential to the flxed horizontal boundaries. 



D. D. Holm k M. F. Staley 



Singular Wave Fronts 137 



[4] The interfacial velocities Ki appear in the variational derivatives of the 
Lagrangian i explicitly as 

Ui + 7vV-f(2/i, + /i,+i) 



PiDi Sui Di 6 

1 D 
+ -(/ii + %i)V/im + v5^-^(/i, + /i,+i), 



1 6e 

Pi SDi 



1 



[5] S. V. Bazdenkov, N.N. Morozov and O.P. Pogutse, 1985, Dispersive effects 
in two-dimensional hydrodynamics. Dokl. Akad. Nauk SSSR 293, 818-822; 
[Sov. Phys. Dokl. 32, 262-264]. 

[6] R. Camassa and D. D. Holm, An Integrable Shallow Water Equation with 
Peaked Solitons. Phys. Rev. Lett. 71 (1993) 1661-1664. 

[7] R. Camassa, D. D. Holm and C. D. Levermore, Long-time effects of bot- 
tom topography in shallow water. Physica D 98 (1996) 258-286. 

[8] J. G. Charney and P. G. Drazin, Propagation of planetary scale distur- 
bances from the lower into the upper atmosphere. J. Geophys. Res. 66 
(1961) 83-109. 

[9] S.Y. Chen, C. Foias, E.J. Olson, E.S. Titi and S. Wynne, The Camassa- 
Holm equations as a closure model for turbulent channel and pipe flows. 
Phys. Rev. Lett. 81 (1998) 5338-5341. 

[10] W. Choi and R. Camassa, Long internal waves of finite amplitude. Phys. 
Rev. Lett. 77 (1996) 1759. 

[11] W. Choi and R. Camassa, Fully nonlinear internal waves in a two-fluid 
system. J. Fluid Meek. 396 (1999) 1-36. 

[12] Desbrun, M., A. N. Hirani, M. Leok and J. E. Marsden [2003], Discrete 
exterior calculus, (preprint). 



D. D. Holm & M. F. Staley 



Singular Wave Fronts 138 



[13] H. R. Dullin, G. A. Gottwald and D. D. Holm, An integrable shallow 
water equation with linear and nonlinear dispersion. Phys. Rev. Lett. 87 
(2001) 194501-04. 

[14] H. R. Dullin, G. A. Gottwald and D. D. Holm, Camassa-Holm, 
Korteweg-de Vries-5 and other asymptotically equivalent equations for 
shallow water waves. Fluid Dyn. Res. 33 (2003) 73-95. 

[15] H. R. Dullin, G. A. Gottwald and D. D. Holm, On asymptotically equiv- 
alent shallow water wave equations. Physica D 190 (2004) 1-14. 

[16] C. Foias, D. D. Holm and E. S. Titi, The Navier-Stokes-alpha model of 
fluid turbulence. Physzca D 152 (2001) 505-519. 

[17] M.F. Gobbi, J.T. Kirby and G. Wei, A fully nonlinear Boussinesq model 
for surface waves. Part 2. Extension to 0{khf. J. Fluid Mech. 405 (2000) 
181-210. 

[18] A. E. Green and P. M. Naghdi, A derivation of equations for wave prop- 
agation in water of variable depth. J. Fluid Mech. 78 (1976) 237-246. 

[19] A. N. Hirani, Discrete Exterior Calculus. Caltech Ph.D. thesis, 2003. 
|http:/7resolver.caltech;edu/CaltechETD:etd-05 202003-095403] 

[20] D. D. Holm, Hamiltonian structure for two-dimensional hydrodynamics 
with nonlinear dispersion. Phys. Fluids 31 (1988) 2371-2373. 

[21] D. D. Holm, Variational multilayer Green-Naghdi fluid equations. In 
preparation. 

[22] D. D. Holm and J. E. Marsden, Momentum Maps and Measure-valued 
Solutions (Peakons, Filaments and Sheets) for the EPDiff Equation. In 
The Breadth of Symplectic and Poisson Geometry, (Marsden, J. E. and 
T. S. Ratiu, eds), Festshrift for Alan Weinstein. Birkhauser Boston, to 
appear. 

[23] D. D. Holm, J. E. Marsden and T. S. Ratiu, The Euler-Poincare equa- 
tions and semidirect products with applications to continuum theories. 
Adv. m Math., 137 (1998) 1-81 

[24] D. D. Holm, J. E. Marsden, T. Ratiu and A. Weinstein, Nonlinear sta- 
bility of fluid and plasma equilibria. Physics Reports 123 (1985) 1-116. 



D. D. Holm & M. F. Staley 



Singular Wave Fronts 139 



[25] D. D. Holm, J. Munn and S. N. Stechmann, Reduced singular solu- 
tions of EPDiff equations on manifolds with symmetry.In preparation for 
Nonlinearity. 

[26] D. D. Holm, V. Putkaradze and S. N. Stechmann, Rotating concentric 
peakons. Submitted to Nonlinearity. 

[27] D. D. Holm, T. Ratnanather, A. Trouve and L. Younes, Soliton dynam- 
ics in computational anatomy, Neurolmage (2004), To appear. 

[28] D. D. Holm and M. F. Staley, Wave Structures and Nonlinear Balances 
in a Family of Evolutionary PDEs. SIAM J. Appl. Dyn. Syst. 2 (3) (2003) 
323-380. 

[29] M.-K. Hsu and A. K. Liu, Nonlinear Internal Waves in the South China 
Sea. Canadian J. Rem. Sens. 26 (2000) 72-81. 

[30] J. K. Hunter and R. H. Saxton, Dynamics of director fields. SIAM J. 
Appl. Math. 51 (1991) 1498-1521. 

[31] J. M. Hyman and M. Shashkov, The Adjoint Operators for the Natural 
Discretizations of the Divergence, Gradient, and Curl on Logically Rect- 
angular Grids, IMACS J. Appl. Num. Math., 25, pp. 413-442 (1997). 

[32] Institute for Mathematics and its Applications (IMA) 
"Hot Topics" Workshop: Compatible Spatial Discretiza- 

tions for Partial Differential Equations, May 11-15, 2004. 
|http : / / www . ima . umn . edu / complex / spring / discretization . html[ 

[33] B. Khesin and G. Misiolek, Euler equations on homogeneous spaces and 
Virasoro orbits. Adv. Math. 176 (2003) 116-144. 

[34] S. Kouranbaeva, The Camassa-Holm equation as a geodesic flow on the 
diffeomorphism group. J. Math, Phys. 40 (1999) 857-868. 

[35] M. Leok, Foundations of Computational Geometric Mechanics. Cal- 
tech Ph.D. thesis, 2004. http://resolver.caltech.edu/CaltechETD:etd-] 
03022004-000251. 



[36] Lew, A., J. E. Marsden, M. Ortiz and M. West [2003], Asynchronous 
variational integrators. Archive for Rat. Mech. An. 167, 85-146. 



D. D. Holm k M. F. Staley 



Singular Wave Fronts 140 



[37] S. T. Li, private commmunication. 

[38] R. Liska, L. Margolin and B. Wendroff, Nonhydrostatic two-layer models 
of incompressible flow. Computers Math. Applies. 29 (1995) 25-37. 

[39] R. Liska and B. Wendrofi^, Analysis and computation with stratified fiuid 
models. J. Comp. Phys. 137 (1997) 212-244. 

[40] A. K. Liu, Y. S. Chang, M.-K. Hsu, N. K. Liang, Evolution of nonlinear 
internal waves in the east and south China Seas. J. Geophys. Res. 103 
(1998) 7995-8008. 

[41] P. J. Lynctt and P. L.-F. Liu, A two-dimensional, depth-integrated 
model for internal wave propagation over variable bathymetry. Wave Mo- 
tion 36 (2002) 221-240. 

[42] P. A. Madsen and Y. Agnon, Accuracy and convergence of velocity for- 
mulations for water waves in the framework of Boussinesq theory. J. Fluid 
Meeh. 477 (2003) 285-319. 

[43] L. G. Margolin, B. T. Nadiga and P. K. Smolarkiewicz, Various approx- 
imations of the shallow flow over an obstacle. Phys. Fluids 1996. 

[44] J. E. Marsden and S. Shkoller Global well-posedness for the Lagrangian 
averaged Navier-Stokes (LANS— a) equations on bounded domains. Phil. 
Trans. R. Soe. Lond. A 359 (2001) 1449-1468. 

[45] J. Miles and R. Salmon, Weakly dispersive nonlinear gravity waves. J. 
Fluid Meeh. 157 (1985) 519-531. 

[46] M. 1. Miller, A. Trouve and L. Younes, On the metrics and Euler- 
Lagrange equations of computational anatomy. Ann. Rev. Biomed. Engin. 
4 (2002) 375-405. 

[47] G. Misiolek, A shallow water equation as a geodesic flow on the BottVi- 
rasoro group. J. Geom. Phys. 24 (1998) 203-208. 

[48] D. Mumford, Pattern theory and vision. In Questions Mathematiques 
En Traitement Du Signal et de L'Image, chapter 3, pages 7?13. Institut 
Henri Poincar, Paris, 1998. 

[49] J. Pedlosky, Geophysieal Fluid dynamics. Springer (1987). 



D. D. Holm k M. F. Staley 



Singular Wave Fronts 141 



[50] C. H. Su and C. S. Gardner, Korteweg-de Vries Equation and Gener- 
alizations. III. Derivation of the Korteweg-dc Vries Equation and Burgers 
Equation. J. Math. Phys. 10 (1969) 536-539. 

[51] G. B. Whitham, Variational methods and applications to water waves. 
Proc. R. Soc. London, Ser. A 299 (1967) 6. 

[52] L. Younes, Computable elastic distances between shapes. SIAM J. Appl. 
Math. 58 (1998) 565-586 

[53] J. Ziegenbein, Spatial observations of short internal waves in the Strait 
of Gibraltar. Deep Sea Res. 17 (1970) 867-875. 



